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Tuesday  Morning  Session 
Wei  come/ Introduction  -  2Lt  O.L.  Warner  (AFGL) 

As  coordinator  of  the  two-day  conference,  Lt.  Warner  opened  the  meeting 
on  behalf  of  the  co-sponsoring  agencies,  the  Air  Force  Geophysics 
Laboratory  (AFGL)  and  the  Defense  Mapping  Agency  (DMA).  She  then  introduced 
Dr.  C.  Martin,  Director  of  the  Advanced  Technology  Division,  DMA. 

Opening  remarks  -  Dr.  C.  Martin 

Dr.  Martin  welcomed  participants  and  attendees.  He  made  reference  to  the 
agenda  noting  his  enthusiasm  over  topics  of  data  processing  and  algorithm 
development,  and  the  breaking  from  topics  of  hardware  development.  He 
expressed  the  urgency  of  the  need  to  develop  the  necessary  means  co 
analyze  the  data.  Software  development  has  a  long  lead  time,  therefore 
Dr.  Martin  urged  the  gravity  community  to  share  insights,  applicable 
software,  anything  that  might  keep  individual  software  developers  from 
"reinventing  the  wheel". 

AFGL/DMA  GGSS  Program  Review  -  Mr.  R.  Boryeson  (AFGL) 

Mr.  Borgeson,  program  manager  of  the  Gravity  Gradiometer  Survey  System, 
presented  a  review  of  those  milestones  that  have  been  accomplished,  and 
those  soon  to  be  accomplished.  He  gave  a  brief  description  of  tha  two 
modes  of  the  system:  the  land  based  and  the  airborne  modes,  including 
mockup  models  and  viewgraph  sketches.  He  stated  the  accuracy  reqjirements 
of  the  system  and  showed  diagrams  and  pictures  of  the  actual  threa-axis 
gravity  gradiometer  that  Bell  Aerospace  Textron  is  under  contract  to 
develop  and  build.  Mr.  Boryeson  stated  that  AFGL  and  DMA  are  in  the 
process  of  determining  which  aircraft  would  best  suit  the  needs  oF  the 
GGSS  program.  The  present  choices  are  the  P3,  the  Convair  58U,  and  the  C13U. 

Questi ons: 

Dr.  S.  Hammer  (University  of  Wisconsin)  -  To  what  accuracy  can  we  estimate 
gravity  using  the  gradiometer?  Is  it  a  prediction  or  demonstrated  performance? 

Mr.  Borgeson  -  .18  arc  seconds,  .9  milligals  and  it  is  a  specification 
and  predicted  accuracy. 

Mr.  M.  May  (NADC)  -  What  navigation  system  would  be  used  on  the  P3,  and 
will  a  gravimeter  be  used? 

Mr.  Boryeson  -  GPS  navigation  and  no  gravimeter. 
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Review  of  1983  Moving  Base  Gravity  Gradiometer  Activities  at  Bell  Aerospace 
Textron  -  Mr.  h»  Metzger  (Bell  Aerospace  Textron) 

Mr.  Metzger  described  the  two  programs  Bell  is  working  on:  the  Navy 
Gravity  Survey  System  (GSS)  and  the  AFGL/DMA  GGSS  Program.  He  described 
the  performance  of  the  Advanced  Development  Model  (ADM)  #1  (Navy)  aboard 
the  DSNS  Vanguard.  Having  been  in  operation  since  February  1982  the 
instrument  failed  less  often  than  expected  and  met  the  feasibility 
requirements  in  all  tests.  The  ADM  #2  has  been  in  operation  since 
September  1983.  Mr.  Metzger  outlined  the  status  and  future  plans  of  ADM 
#2  stating  that  tne  objective  was  to  provide  a  test  bed  for  verifying 
modifications  considered  for  operational  systems,  iince  Mr.  Borgeson  had 
already  covered  the  physical  or  hardware  portions  of  the  GGSS,  Mr.  Metzger 
focused  his  attention  on  a  discussion  of  error  analysis  and  track  spacing. 

Questions: 

Dr.  H.J.  Paik  (University  of  Maryland)  -  What  is  the  origin  of  the  red  noise? 

Mr.  Metzger  -  Drift  in  the  even  order  calibration  coefficients  of  the 
accelerometers. 

Dr.  Paik  -  Is  the  noise  proportional  to  1/f? 

Mr.  Metzger  -  It  is  proportional  tr  o(  tween  1/f  and  (l/f)^*^. 


Comparison  of  At-Sea  Gradiometer  Test  Results  With  An  Independent  Gravity 
Gradient  Reference  -  Mr.  D.  Benson,  Jr.  and  Mr.  A.  Zorn  (Dynamics  Research 
Corporation) 

Mr.  B.  Regenauer,  Mr.  Benson,  and  Mr.  Zorn  are  co-authors  of  this  paper 
and  it  was  presented  by  the  latter  two.  Mr.  Benson  discussed  the  generation 
of  a  gravity  gradient  reference  off  the  north  coast  of  Puerto  Rico  and 
subsequent  comparison  of  Gravity  Gradiometer  gradients  with  the  reference. 

The  gradient  reference  was  derived  using  position  differences  between  an 
accurate  marine  inertial  system  and  an  accurate  geodetic  position  reference. 
Autotape.  Similar  techniques  have  been  used  for  vertical  deflections 
over  land  but,  in  contrast,  require  frequent  stopping  for  gyro  recalibration. 
Mr.  Zorn  presented  part  two,  the  feasibility  of  the  at-sea  survey  technique 
which  was  established  using  covariance  error  analysis.  Subsequently  a 
test  plan  was  developed  to  acquire  data  for  map  development,  and  the 
survey  was  performed  with  equipment  onboard  the  DSNS  Vanguard.  Returns 
over  the  same  paths  and  crossing  paths  during  the  survey  allow  separation 
of  time  varying  gyro  and  acclerometer  errors  from  position  varying 
gradients  and  deflections  when  the  survey  data  is  processed.  The  analysis 
and  data  processing  include  local  geodetic  constraints  on  the  gradients 
to  satisfy  the  conditions  of  the  earth's  gravitational  potential.  Mr. 

Zorn  concluded  by  mentioning  that  the  survey  paths  were  repeated  on  a 
later  at-sea  test  with  the  Bell  Gravity  Gradiometer  onboard.  Gravity 
gradients  from  the  gradiometer  compared  remarkably  well  with  the  gravity 
gradient  reference  along  the  survey  paths. 

Questions: 

Mr.  Borgeson  -  Is  the  program  under  contract? 

Mr.  Benson  -  Yes,  under  contract  with  Sperry  Management  Systems,  SP24. 
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Mr.  S.  Jordan  (Geospace)  -  Can  you  give  an  actual  accuracy  for  test  results? 
Mr.  Benson  -  No,  it  is  classified. 

Mr.  P.  Fell  (DMAHTC)  -  What  is  the  Autotape  system? 

Mr.  Benson  -  A  radar  navigation  line  of  sight  positioning  system. 

Mr.  W.  Gumert  (Carson  Geoscience)  -  What  speed  were  you  survey! ni  at? 

Mr.  Zorn  -  5  to  8  knots. 

Dr.  W.  Heller  (TASC)  -  Is  f^ere  an  RMS  difference  between  map  ano  measured 
gradients? 

Mr.  Benson  -  They  exist  but  are  classified. 

A  discussion  followed  between  Dr.  Heller  and  Mr.  Benson  concerning 
possible  RMS  differences  and  how  they  depend  on  data  processing  ind  use  of 
low  frequency  information. 


Current  GGSS  Performance  Expectation  and  Error  Analysis  Allocation  - 
Dr.  W.  Heeler  (TASC ) 

Dr.  Heller  opened  with  a  review  of  test  and  survey  geometry  of  tie  GGSS 
and  the  desired  goals.  He  discussed  the  error  sources  as  well  a;  the 
GGSS  sensitivity  to  variation  in  the  frequency  content  of  the  grivity 
field.  Mentioned  as  well  were  the  contribution  of  each  class  of  errors 
and  the  sensitivity  to  survey  altitude.  In  his  conclusions  Dr.  Heller 
stated  that  the  bottom-line  is  that  the  performance  of  the  GGSS  is  expected 
to  fall  well  within  prescribed  specifications. 

Quest! ons; 

Mr.  J.  Brozena  (NRL)  -  The  downward  continuation  and  truncation  errors 
refer  to  the  center  of  the  survey.  What  happens  as  you  go  away  From  the 
center? 

Dr.  Heller  -  We  have  to  be  concerned  with  edge  effects.  Testing  will  be 
done  at  the  center  of  the  survey  region  to  minimize  edge  effects. 

Mr.  Zorn  -  Was  only  one  gradient,  Tzz,  used? 

Dr.  Heller  -  Yes. 

Mr.  Zorn  -  What  if  other  gradients  are  used,  such  as  those  that  are  more 
directly  related  to  the  vertical  deflection  by  integration?  The  error 
could  be  overstated  because  these  gradients  are  not  used. 


Dr.  Heller  -  We  used  Breakwell's  results  that  the  best  yradient  to  use  is 
Tzz.  We  incurred  a  penalty  of  a  factor  of  2  by  not  using  other  gradients 
and  the  errors  in  Tzz  were  correspondingly  scaled. 

Mr.  M.  Molny  (Speiry)  -  Why  not  put  in  a  wider  grid  spacing  for  the 
outside  areas? 

Dr.  Heller  -  That  s  a  good  introduction  to  the  next  talk! 

Dr.  Jordan  -  Is  tie  downward  continuation  to  the  physical  surface  or  to  a 

reference  surface.’ 

Dr.  Heller  -  This  has  to  be  resolved.  We  don't  want  to  continue  the  data 
below  the  surface.  For  convenience  we  chose  the  physical  surface. 

Or.  Jordan  -  I  doi't  iielieve  we  can  meet  the  goals  of  the  program  using 

the  physical  surf  ice  is  a  reference.  We  can't  downward  continue  in  a 
mountainous  area. 

Benson,  Jordan,  aid  H-'ller  further  dicussed  the  choice  of  the  reference 
surface. 

Dr.  H.  Baussus  vo  i  Luetzow  (ETL)  -  We  .leed  to  consider  the  shortcomings 
of  least  squares  .ollocation  at  high  frequencies  (the  use  of  stationary 
models  could  intrjduce  errors  of  1U-2U%).  No  terrain  effects  have  to  be 
considered. 

Dr.  Heller  -  We  plan  to  survey  over  smooth  areas  to  test  and  demonstrate 
the  gradiometer  and  should  be  able  to  handle  mountainous  terrain  in  the 
future. 

Mr.  M.  Trageser  -  Why  do  we  need  to  downward  continue  at  all?  Why  not 
estimate  gravity  at  altitude? 

Or.  Heller  -  We  want  point  estimates  at  the  surface. 

There  was  further  discussion  on  this  question  between  Heller,  Trageser, 
and  Martin. 


Tuesday  Afternoon  Session 


Efficient  Proces'.ing  of  Gradiometer  Output  Using  a  Kalman  Filter  - 
Dr.  J.  Hutcheson  (Bell  Aerospace  Textron) 

Dr.  Hutcheson  male  several  points  concerning  the  processing  of  the  data 
to  be  collected  oy  the  gravity  gradiometer  system.  First,  there  are  just 
too  many  data  points;  they  must  be  cort^ressed.  He  presented  a  method  to 
process  gradient  data  using  a  straight  line  Kalman  filter.  Dr.  Hutcheson 
showed  comparisons  of  the  STAG  model  and  state  space  approximations.  He 
compared  results  of  fixed  point  smoothing  and  least  squares  collocation. 
Results  were  simulated  using  a  MASCON  synthetic  gravity  field. 
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Questions: 

Unidentified  questioner  -  What  wavelengths  were  modeled  with  the  MASCON 
model ? 

Or.  Hutcheson  -  The  North  Atlantic  Gravity  Model  was  used.  The  v^avelengths 
are  possibly  classified. 

Ur.  C.  dekeli  (AFGL)  -  When  using  the  straight  line  filter,  are  correlations 
from  track  to  track  taken  into  account? 

Dr.  Hutcheson  -  Two  dimensional  correlations  are  taken  into  accoi  nt  in  a 
subsequent  stage  of  processing. 

Mr.  P.  Ugincius  (NSWC)  -  What  about  estimating  away  from  the  track? 

Dr.  Hutcheson  -  One  can  always  contract  the  estimation  point  to  the 
nearest  grid  point. 


On-going  Work  in  Post  Mission  Data  Processing  Algorithm  Development  - 

DF.  ¥.“nener~rn\'?n) -  - — - — 

Ur.  Heller  presented  some  background  information  on  the  GGSS  program  and 
described  what  will  be  required  for  post  mission  analysis  of  the  data 
collected.  He  discussed  various  survey  and  data  processing  issues  as 
well  as  prefilter  design  considerations.  On  this  topic  Dr.  Heller  concludes 
that  a  12  second  filter  time  constant  will  be  needed.  On  the  toeic  of 
estimation  approacnes  to  the  gravity  disturbance  vector.  Dr.  Heller 
proposed  two  alternatives: 

a)  a  sequential  complementary  filter. 

b)  an  optimized  template  algorithm. 

Numerical  data  showing  applications  of  these  alternatives  were  presented. 

Dr.  Heller  summarized  by  stating  that  on-going  development  is  not  a  neat 
"packaged  deal"  but  an  active  process.  He  then  proceeded  to  outline  the 
development  of  the  algorithms. 

Questi ons : 

Mr.  Zorn  -  How  sensitive  is  the  estimation  if  the  aircraft  is  net  able 
to  follow  the  regular  grid? 

Dr.  Heller  -  The  system  noise  includes  navigation  error.  If  the  navigation 
error  becomes  too  great,  we  may  cease  data  collection. 

Dr.  b.  Hammer  -  How  did  you  simulate  the  test  data? 

Or.  Heller  -  With  ocean  trench  data  by  computing  gradients  which  were 
reinverted  creating  an  anomaly  field. 

Mr.  Molny  -  What  is  the  actual  mapping  area? 

Dr.  Heller  -  320  X  32U  kilometers. 

Mr.  Horgeson  -  There  are  two  test  sites.  Central  Florida  and  Cheyenne, 
Wyoming,  are  being  considered  based  on  weather  data  and  availabe  gravity 


Mr.  Molny  -  When  cor'densiny  lata,  averayiny,  etc.,  do  you  consider  the 
survey  design? 

L)r.  Heller  -  We  must  make  sure  the  data  are  regular. 

Mr.  Molny  -  One  could  enlarge  the  size  of  the  square  by  having  laryer 
track  intervals  in  the  outside  areas  to  pick  up  larger  wavelength 
i nformati on. 

Ur.  Heller  -  We  may  do  th,s  ror  the  test  if  we  are  interested  only  in 
gravity  estimates  in  the  inner  zone. 

Uyincius,  Hutcheson,  Fell,  Zorn,  Benson,  and  Heller  discussed  use  of  the 
template  method. 


Three-axis  Super  ronciuctiny  Gravity  Gradiometer  -  Or.  H.J.  l-’aik  (University 
of  Maryland) 

Ur.  Paik  opened  his  presentation  with  a  review  cf  his  program  to  develop 
a  single  axis  version  of  a  cryogenic  gravi’y  gradiometer.  Dr.  Paik  then 
went  on  to  describe  plans  to  develop  a  th^ee-axis  superconducting  gravity 
gradiometer.  He  elaborated  mathematically  on  the  expected  sensitivity  of 
the  instrument  and  explained  its  physical  design  with  slides  and  viewgraphs. 
Or.  Paik  described  i.he  new  pendulum  suspension  and  error  compensation. 

For  adequate  error  compensation  Dr.  Paik  proposed  the  development  of  a 
six-axis  superconducting  accelerometer  and  described  in  detail  its 
configuration,  design  features,  and  sensitivity  as  well  as  its  operation 
with  a  gradiometer.  He  concluded  with  a  proposed  time  line  for  the 
development  of  both  a  three-axis  GGM  and  a  six-axis  accelerometer.  His 
ultimate  aspirations  include  space  born  testing  of  the  GGM  in  an  ea^th 
orbit  aboard  the  space  shuttle  by  the  year  1990. 

Questions: 

Ur.  Heller  -  At  what  time  will  an  eror  of  1U~^  E/(riz)^/^  be  demonstrated? 

Ur.  Paik  -  This  is  difficult  to  predict.  We  may  have  to  do  temperature 
control  experiments  to  yet  rid  of  squid  noise.  We  had  hoped  to  have  had 
data  2  years  ago,  but  we  need  to  improve  the  platform.  With  a  little 
confidence  and  if  the  platform  is  UK,  maybe  next  week.  We  are  trying  to 
yet  better  noise  levels  in  the  platform.  Maybe  next  year. 

Mr.  Hastings  (Sperry)  -  If  you  were  to  add  proportional  feed  back,  would 
it  improve  the  accuracy? 

Ur.  Paik  -  This  would  not  improve  the  accuracy.  It  would  increase  the 
dynamic  range. 

Dr.  Jordan  -  GHM  will  get  100  km  resolution  at  a  cost  of  1/3  billion 
dollars.  Would  you  care  to  comment  on  the  resolution  and  cost  of  the 
satellite  gradiometer? 

Or.  Paik  -  For  lir^  E  sensitivity  we  get  lUU  km  resolution.  To  improve 
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the  resolution  by  a  factor  of  two  we  need  to  improve  the  gradiomere''  by  a 
factor  of  lUU. 


Overview  of  Airborne  Gravity  and  Comparison  of  Gravity  and  Gradie nt 
AnomalTes  -  Dr.  S.  Hammer  (University  of  Wisconsin) 

Dr.  Hammer  introduced  his  topic  with  a  very  brief  history  of  gravity 
measurements  in  a  moving  vehicle  and  a  description  of  the  Carson  Helicopter 
Gravity  Measuring  System  (HGMS).  HGMS  has  been  operational  for  about  four 
years  and  observes  an  average  of  SC  JO  km  of  line  data  per  month  vith  a 
probable  error  of  the  order  of  1  mgal  and  anomaly  resolution  limit  of 
about  3  km.  Survey  altitude,  as  specified,  can  vary  from  low  terrain 
clearance  up  to  about  14,000  ft  (4300  m)  above  sea  level  and  over  any 
type  of  terrain  including  environmentally  restricted  areas  such  cs  the 
Wildlife  Reserve  in  Alaska.  Ir,  his  presentation  Or.  Hammer  compr  red  the 
magnitudes  and  rates  of  attenuation  with  altitude  of  gradient  ant  gravity 
anomalies.  He  emphasized  that  the  rate  of  attenuation  of  all  components 
of  the  gravity  gradient  is  the  same  for  a  given  type  of  anomaly.  The 
gradient  attenuates  .,'*h  altitude  an  order  of  magnitude  faster  than  the 
corresponding  gravity  anomaly.  In  conclusion  Or.  Hammer  raised  'he 
question:  to  accomplish  the  basic  objectives  of  this  conference,  what  is 
the  best  way  to  proceed?  (a)  Concentrate  efforts  on  the  development  of 
theory  and  instrumentation  for  gravity  gradiometry;  or  (b)  Undertake  a 
world-wide  gravity  mapping  program  with  the  rapid  and  precise  HGMS  whi  ih 
i s  now  avai 1 abl e. 

Questions:  (deferred  to  Mr.  W.  Gumert  of  Carson  Geoscience) 

Mr.  Borgeson  -  Did  you  have  any  problems  with  civil  air  traffic? 

Mr.  Gumert  -  There  were  no  problems,  as  long  as  we  filed  a  fligit  plan 
and  flew  low. 

ijr.  Heller  -  How  far  out  to  sea  die.  you  fly? 

Mr.  Gumert  -  150-180  km  using  line  of  sight  radar. 

Mr.  Molny  -  How  do  you  separate  vertical  acceleration  from  gravity  anomalies? 
Ur.  Hammer  -  Using  barometric  and  radar  altimeters. 


Wednesday  Horning  Session 

(Due  to  unexpected  circumstances,  the  order  of  presentation  as  shown  on 
the  agenda  was  slightly  rearranged.) 

A  Simpl i f i cati or  of  the  Least  Squares  Determination  of  the  Gravity  Field 
from  Low-low  SaLellite  Tracking  Data  Utilizing  An  Operator  Representation 
of  the  Range  Observable  -  Dr.  V.  Reinhardt  (Bendix  Field  Enginef ring 
Corporation ) 

Dr.  Reinhardt  introduced  his  paper  by  explaining  that  one  of  the  purposes 
of  NASA/Goddard  Space  Flight  Center's  proposed  Geopof^ntial  Research 
Mission  is  to  map  the  gravity  field  of  the  Earth  to  1-11  mgal  in  41,Eb3 
1  degree  by  1  degree  blocks  using  range  rate  tracking  data  between  two 


polar  low  orbitiny  satellites  in  the  same  orbit.  Because  there  are 
41,2b3  blocks  in  the  yravity  map,  it  takes,  in  yeneral,  at  least  9,4x10^2 
real  floating  point  operations  (flop)  to  solve  the  least  squares  fit 
(LSQF)  matrix  equation  used  to  determine  the  gravity  potential  from  the 
tracking  data.  Dr.  Reinhardt  showed  that  with  the  circular  orbit  infinite 
trajectory  length  approximation,  and  an  iterative  method,  one  obtains  an 
approximate  covariance  matrix  for  the  coefficients  in  less  than  l.BxlO^'^ 
flop.  In  his  presentation  Dr.  Reinhardt  also  demonstrated  the  application 
of  this  method  to  airborne  gradiometry. 


WGS  Gravity  Gradient  Reference  Ellipsoid  and  Earth  Gravitational  Model  - 
Mr.  M.  May  (NADC) 

Mr.  May  described  the  gravity  gradient  ellipsoid  in  terms  of  geographic 
latitude  using  WGS-7i?  reference  ellipsoid  parameters.  Additionally,  he 
described  a  representation  of  earth's  gravity  gradients  by  a  Fourier 
series  for  wavelengths  down  to  120  nautical  miles. 

Questions : 

Dr.  Hammer  -  What  wavelength  ranges  were  used? 

Mr.  May  -  IHU  X  130  degree  field  whicii  cor'-esponds  to  wavelengths  df'wn  to 
120  nautical  miles. 

Dr.  Heller  -  Comment:  for  stochastic  modeling,  models  have  been  determined 
for  the  degree  variances  of  the  gravity  field. 

Dr.  Jekeii  -  What  is  the  advantage  of  using  gradient  maps  over  gravity  maps 

Mr.  May  -  It  is  advantageous  to  mix  data  for  vertical  deflection 
determination  at  the  instrument  (gradient)  level. 


Re-looking  at  MASCONS  for  Representing  Gravity  Survey  Data  - 
Or.  J.  Breakwel'l  (Stanford  University) 

Or.  Breakwell  introduced  his  presentation  with  the  statement  that  the 
representation  of  the  higher  harmonic  part  of  the  Earth's  gravity  field 
by  a  finite  grid  of  MASCONS  in  two  layers,  one  at  the  Earth's  surface  and 
one  at  the  bottom  of  the  Earth's  crust,  leads  to  an  error  in  the  gravity 
components  at  altitude  (in  addition  to  the  error  from  the  estimation  of 
density  fluctuations  in  the  two  layers).  His  discussion  examined  these 
errors  for  various  ratios  of  grid  size  to  gradiometer  altitude. 


Integral  Formulas  Relating  Gravity  Gradients  to  Vertical  Deflections  - 
Mr.  A.  Zorn  (Dynamics  Research  Corporation) 

Mr.  Zorn  presented  a  number  of  integral  formulas  relating  gravity 
gradients  to  verticdl  deflections.  The  formulas  are  derived  from  standard 
geodesy  theory  and  are  similar  to  the  Stokes  and  Veni ng-Mei nesz  integral 
formulas.  He  explaineu  that  the  deterministic  formulas  provide  a  physical 
explanation  of  certain  properties  common  to  various  self-consistent 
stochastic  gravity  models.  For  example,  the  high  degree  of  correlation 
between  the  gravity  gradient,  Txz,  and  the  vertical  deflection  Tx,  in 
most  gravity  models  can  be  physically  justified  based  on  a  deterministic 
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inteyral  formula.  The  integral  formulas  may  lead  to  insight  into  the 
construction  of  gradiometer  filters  to  estimate  vertical  deflections 
which  are  insensitive  to  stochastic  gravity  models.  Mr.  Zorn  discussed 
implementation  considerations  involving  discretization  and  truncation  of 
the  integral  formulas  using  surface  gravity  data  over  the  Blake  Escarpment. 
He  concluded  with  several  statements  including  that  vertical  deflection 
maps  can  be  const''ucted  from  gradiometer  survey  data  using  Stokes'  formula. 

Questions; 

Mr.  Ugincius  -  Going  from  Tzz  one  can  get  any  other  gradient  quantity. 

Why  measure  these  other  components? 

Mr.  Zorn  -  To  get  other  components  we  need  Tzz  over  the  entire  earth's 
surface.  We  therefore  measure  the  other  components  to  get  more  nformation 
especially  in  limited  areas. 

Mr.  Ugincius  -  Would  the  tree  diagram  hold  for  spherical  formula  :ions? 

Mr.  Zorn  -  Yes 

Ur.  D.  UeBra  (Stanford  University)  -  Are  some  paths  more  sensitize? 

Mr.  Zorn  indicated  on  his  diagram  certain  paths  that  would  be  affected  by 
low  frequency  errors. 

Mr.  May  -  The  tree  diagram  can  be  derived  from  spherical  harmonic 
expansion.  Is  there  distinction  between  gravity  anomaly  and  gravity 
di sturbance? 

Mr.  Zorn  -  The  difference  is  negligible. 

Mr.  A.  Rufty  (NSWC)  -  If  we  have  gradient  information  only,  we  need 
integration  constants  to  get  gravity  quantities. 

Zorn  -  We  would  need  low  frequency  information  to  obtain  the  constant  of 
integration. 


Real-Time  Vertical  [Reflection  Estimation  Using  Gradient  Data  - 
Ur.  W.  Feldman  (Sperry;  Co-author,  Mr.  B.  Epstein) 

Ur.  Feldman  discussed  techniques  and  algorithms  for  making  real-time 
estimates  of  vertical  deflections  from  gravity  gradient  measurements. 

First  he  outlined  the  techniques  and  algorithms  used  to  condition  sensor 

data  for  effects  of  high-frequency  noise,  pressure  sensitivity,  self- 

gradients,  bias/trend,  and  carouselling.  He  then  showed  plots  cif  typical  ^ 

gradient  data  before  and  after  preprocessing  to  illustrate  the  ( ffecti veness 

of  the  preprocessing.  Second,  a  batch  filter  formulation  for  estimating 

deflections  from  gradient  and  SEASAT  map  data  was  illustrated.  Finally, 

Dr.  Feldman  demonstrated  filter  results  by  plots  of  filter  estimates  of 
vertical  deflection  vs.  reference  values  for  several  ship  tracks. 

Questions: 

Dr.  Heller  -  Why  pull  out  the  reference  field  in  the  instrument  frame 
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rather  than  in  the  NED  frame? 

Dr.  Feldman  -  It  is  equivalent. 

Dr.  Heller  -  What  are  the  practical  considerations? 

Dr.  Feldman  -  It  was  simpler  to  do  in  the  instrument  frame  and  there  was 
less  chance  of  error. 


Gravity  Gradiometer  Application  to  Tunnel  Detection  -  Mr.  Jircitano  (Bell) 
Mr.  Jircitano  described  a  project  being  conducted  by  the  Arrriy  in  the 
field  of  tunnel  detection  by  means  of  measuring  differences  in  gravity 
gradients.  There  are  two  tunnel  detection  problems:  1)  the  detection  of 
new  tunnels  (those  which  had  been  dug  after  a  previous  survey  of  the 
area)  2)  the  detection  of  old  tunels  (those  that  were  present  before 
any  previous  survey). 

Quest i ons: 

Dr.  Reinhardt  -  Whac  is  the  feasibility  of  tunnel  detection  if  they  were 
masked  by  high  density  fill? 

Mr.  Jircitano  -  The  possibility  exists  but  to  compensate  one  would  have 
to  half  fill  the  tunnel  with  steel  which  is  impractical. 

Dr.  DeBra  -  Seasonal  groundwater  has  a  significant  mass.  Is  this  a 
problem  for  new  tunnel  detection? 

Mr.  Jircitano  -  Could  be,  but  water  would  not  collect  in  a  tunnel  formation. 


Closing  Remarks  -  Lt.  D.  Warner 

Lt.  Warner  thanked  everyone  for  attending  this  year's  conference,  especially 
those  who  presented  papers.  She  also  offered  a  special  thanks  to  her  AFGL 
colleagues  who  helped  with  audio/visual  equipment  and  note  taking.  As  a 
date  for  next  year's  conference  had  not  yet  been  established,  she  assured 
all  that  they  world  be  informed  of  the  date  as  soon  as  it  was  reserved. 
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SUMMARY  ERROR  ANALYSIS 
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GGSS  FUNCTIONAL  BLOCK  DIAGRAM 
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LAND  VEHICLE  INSTALLATION 


Bell  Aerospace 


66SS  BASE  STATION  SUPPORT  EQUIPMENT 
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PLATFORM  &  ENCLOSURE  -  TOP  VIEW 


PLATFORM  &  ENCLOSURE  -  SIDE  VIEW 


Minimal  Changes  From  ADM  Platform 
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Terminal  Block  of  ADM  Azimuth  Gimbal  Replaced  by  ^  Connectors  as  on  FSD 
All  Major  Components  are  Drawn  and  Being  Checked 


AIRCRAFT  ACCELERATION  MEASUREMENTS 
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A  Shock  Level  Greater  Than  20  gs. 

Cabinets  Will  be  in  EQUIPTO  Racks  -  Two  Sizes. 


Completed  Design  for  GGSS  and  GSS  Provides  Flexibility  for  Slip  Ring 

Selection  and  Model  VII  Pendulosity. 


GGS  UNINTERRUPTABLE  POWER  SUPPLY  (UPS) 
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SCHEMATIC  ILLUSTRATION  OF  ROTATING  FIXTURE 
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GGI  SIGNAL  PROCESS  BLOCK  DIAGRAM 
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BOBBIN  AND  LEG  ASSEMBLY  -  LIGHT  PROOF 


ACCELEROMETER  ERROR  MODEL 
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GGI  OUTPUT  ERRORS  AS  A  FUNCTION  OF  ACCELERATION 
DISTURBANCES  ERROR  COEFFICIENTS  AND  ATTITUDE 
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PLATFORM  PERTURBATION  WITH  AXIAL  SHAKE  LOOP  SETTING 
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COMPARISON  OF  AT-SEA  GRADIO METER  TEST  RESULTS 
WITH  AN  INDEPENDENT  GRAVITY  GRADIENT  REFERENCE 


Donald  O.  Benson,  Jr. .  Alan  H.  Zorn,  and  Bernard  J.  Regenauer 
Dynamics  Research  Corporation 
60  Concord  Street 
Wilmington,  Massachusetts  01887 

! 

At-sea  demonstration  and  testing  of  the  Bell  Gravity  Gradiometer  near 
coastal  regions  is  desirable  to  allow  access  of  engineering  personnel  for 
equipment  grooming  and  maintenance,  and  to  obtain  reasonable  signal -to - 
noise  ratios  since  large  deflections  and  gradients  occur  near  many  coastal 
regions  (Ref.  1).  Reference  gravity  maps  in  these  regions  either  do  not 
exist,  do  not  have  sufficient  accuracy  when  derived  from  satellite  altimetry 
or  are  extremely  time  consuming  to  survey  using  gravimetric  methods. 

This  presentation  discusses  the  generation  of  a  gravity  gradient 
reference  off  the  north  coast  of  Puerto  Rico  and  subsequent  comparison  of 
Gravity  Gradiometer  gradients  with  the  reference.  The  gradient  reference 
was  derived  using  position  differences  between  an  accurate  marine  inertial 
system  and  an  accurate  geodetic  position  reference.  Autotape.  Similar 
techniques  have  been  used  for  vertical  deflections  over  land  but,  in  contrast, 
require  frequent  stopping  for  gyro  recalibration  (Ref.  2).  Feasibility  of  the 
at-sea  survey  technique  was  established  using  covariance  error  analysis. 
Subsequently  a  test  plan  was  developed  to  acquire  data  for  map  development, 
and  the  survey  was  performed  with  equipment  onboard  the  USNS  Vanguard. 

Returns  over  the  same  paths  and  crossing  paths  during  the  survey  allow 
separation  of  time  varying  gyro  and  accelerometer  errors  from  position 
varying  gradients  and  deflections  when  the  survey  data  is  processed.  The 
analysis  and  data  processing  include  local  geodetic  constraints  on  the  gradients 
to  satisfy  conditions  of  the  earth's  gravitational  potential  (Ref,  3). 

The  survey  paths  were  repeated  on  a  later  at-sea  test  with  the  Bell 
Gravity  Gradiometer  onboard.  Gravity  gradients  from  the  gradiometer 
compare  remarkably  well  with  the  independent  gravity  gradient  reference 
along  the  survey  paths. 
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OVERVIEW  OF  AIRBORNE  GRAVITY  AND  COMPARISON  OF  GRAVITY  AND  GRADIENT  ANOMALIES* 

Sigmund  Hammer 

University  of  Wisconsin,  Madison,  WI  53706 
William  R.  Gumert 

Carson  Geoscience  Co.,  Perkasie,  PA  18944 


ABSTRACT 

The  topic  is  introduced  by  a  very  brief  history  of  gravity  measurements 
in  a  moving  vehicle.  The  Carson  Helicopter  Gravity  Measuring  System  (HGMS) 
is  then  described.  HGMS  has  been  operational  about  four  years  and  observes 
an  average  of  5000  km  of  line  data  per  month  with  a  probable  error  of  the 
order  of  1  mgal  and  anomaly  resolution  limit  of  about  3  km.  Survey  altitude, 
as  specified,  can  vary  from  low  terrain  clearance  up  to  about  14,000  ft 
(4300  m)  above  sea  level  and  over  any  type  of  terrain  including  environmental¬ 
ly  restri.'.  cd  areas  such  as  the  Wildlife  Reserve  in  Alaska. 

The  magnitudes  and  rates  of  attenuation  with  altitude  of  gradient  and 
gravity  anomalies  are  then  compared.  It  is  pertinent  to  emphasize  that  the 
rate  of  attenuation  of  all  components  of  the  gravity  gradient  is  the  same  for 
a  given  type  of  anomaly.  The  gradient  attenuates  with  altitude  an  order  of 
magnitude  faster  than  the  corresponding  gravity  anomaly. 

As  a  conclusion  the  question  is  raised:  to  accomplish  the  basic  objec¬ 
tives  of  this  Conference,  what  is  the  best  way  to  proceed?  (a)  Concentrate 
efforts  on  the  development  of  theory  and  instrumentation  for  gravity  gradi- 
ometry;  or  (b)  Undertake  a  world-wide  gravity  mapping  program  with  the  rapid 
and  precise  HGMS  which  is  now  available. 


♦Prepared  for  presentation  at:  Twelfth  Moving  Base  Gravity  Gradiometry  Re¬ 
view,  U.S.  Air  Force  Academy,  Colorado  Springs,  CO,  February  14-15,  1984 
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OVERVIEW  QF  AIRBORNE  GRAVITY  AND  COMPARISON  OF  GRAVITY  AND  GRADIENT  ANOMALIES 

Sigmund  Hammer  and  William  R.  Gumert 


INTRODUCTION 

So  far  as  I  know  there  are,  as  yet,  no  published  field  data  from  func¬ 
tional  gravity  gradiometers .  Oar  overview,  therefore,  will  assess  the  pres¬ 
ent  state  of  gravity  observations  leading  up  to  the  present  Helicopter  Gravity 
Measuring  System  (HGMS).  This  gravity  overview  may  provide  some  useful  back¬ 
ground  for  evaluating  the  comparative  potential  usefulness  of  airborne  gravity 
versus  airborne  gradiometer  measurements  for  this  symposium. 

First  a  very  brief  history  of  gravity  measurements  in  a  moving  vehicle. 
The  first  measurements  of  gravity  in  a  moving  vehicle  were  accomplished  55 
years  ago  by  Vening-Mei nesz  in  a  submerged  submarine  using  special  three- 
pcnduljr  apparatus.  His  spectacular  results  over  the  Indonesian  Islands  Arc- 
Trench  system  were  early  precursors  of  the  subduction  phenomenon  which  is  a 
major  element  in  today's  plate  tectonics  theory.  Shipborne  gravity  surveys 
with  modern  high-precision  gravimeters  began  40  years  ago  and  continue  today 
on  a  commercial  scale. 

Airborne  gravity  measurements  have  been  under  investigation  since  about 
1958.  Tests  in  fixed-wing  aircraft  have  yielded  enough  encouragement  to  jus¬ 
tify  continuation  but  satisfactory  results  have  not  been  achieved  up  to  now. 
Successful  commercial  gravity  surveys,  using  specially  stabilized  helicop¬ 
ters,  have  been  operational  since  1979.  High  frequency  anomalies,  applicable 
to  petroleum  exploration,  are  being  mapped  at  requested  altitudes  from  low 
terrain  clearance  to  more  than  14,000  feet  (4300  m).  Since  HGMS  is  quite  new 
it  may  be  useful  to  present  a  very  brief  explanation  of  the  operation. 


HELICOPTER  GRAVITY  MEASURING  SYSTEM  (HGMS) 


The  helicopter  measoreraent s  are  made  in  a  large  especially  stabilized 
Sikorski  S-61 .  The  first  slide  shows  the  helicopter  in  flight.  Gravity  is 

SLIDE  #I :  Carson  helicopter  in  flight 

measured  by  a  LaCoste-Roraberg  sea-air  meter  on  a  greatly  improved  stable  plat¬ 
form  at  the  center  of  gravity  in  the  helicopter.  It  may  be  of  interest  to 
point  out  that  the  helicopter  has  been  modified  so  that  the  center  of  gravity 
does  not  shift  as  thousands  of  pounds  of  fuel  are  consumed.  Total  field  mag¬ 
netic  field  data  are  recorded  simultaneously  by  a  trailed  proton  precession 
magnetometer.  Terrain  clearance  is  recorded  by  narrow  bean  laser.  Ail  data 
are  recorded  digitally  each  second.  Flying  is  done  at  night  to  minimize  air 
turbulence.  Flight  speed  is  normally  50  knots  but  can  be  specified  for  the 
purpose  at  hand  up  to  100  knots  (185  km/hr). 

The  procedure  which  brought  success  was  to  fly  a  more  or  less  square 
grid  of  intersecting  lines  as  illustrated  on  slide  2.  In  this  project  the 

SLIDE  yA2:  Flight  pattern 

line  spacing  was  3,300  feet  (1  km)  and  flight  altitude  1000  feet  (305  m).  The 
area  is  a  marshy  jungle  of  extremely  difficult  access  on  the  ground.  Note 
that  north-south  flights  are  avoided  to  improve  the  EHtvHs  corrections. 


THE  ATTENUATION  PROBLEM 


The  principal  elements  in  the  operation  are  well  illustrated  in  a  demon¬ 
stration  survey  over  two  well-known  salt  domes  in  the  Texas  Gulf  Coast  (Car- 

SLIDE  #3:  Gulf  Coast  survey 

son,  1981).  The  survey  consisted  of  20  lines,  spaced  one-half  statute  mile 
(0.805  km \  within  the  area  bordered  by  dashed  lines  as  shown  on  Figure  3. 

The  contoured  results  clearly  define  the  airborne  gravity  ancxnalies  on  the 
two  salt  domes  in  the  shaded  areas.  The  contour  interval  is  1/2  mgal. 

A  comparison  with  ground  data  is  illustrated  on  Figure  4.  The  salt  domes 

SLIDE  PU:  Profile  results 

are  sketched  at  the  bottom  of  the  figure.  Note  the  exaggerated  vertical 
scal'^  on  the  domes  -  the  caprock  is  very  shallow.  Conventional  ground  grav¬ 
ity  data  along  the  location  of  the  flight  line  across  the  centers  of  the  dome 
are  shown  by  the  central  profile.  These  data  were  kindly  provided  by  inter¬ 
ested  oil  companies,  but  not  until  after  completion  of  the  HGMS  project.  The 
gravity  profile  measured  in  flight  is  the  solid  curve  at  the  top  of  the  fig¬ 
ure.  Upward  continuation  of  the  ground  data,  without  filtering,  are  shown  by 

the  dotted  curve  ( . ).  For  comparison  with  the  helicopter  measurements 

the  upward  continued  ground  data  were  then  smoothed  by  the  60-second  filter 
which  removes  high  frequency  motional  acceleration  effects  in  the  HGMS  obser¬ 
vations  ( - ).  The  differences  between  the  airborne  and  filtered  ground 

data  nowhere  exceed  i  mgal.  However,  it  is  obvious  that  the  sharp,  positive 
caprock  anomalies  are  strongly  attenuated.  This  will  be  examined  below. 

A  simulated,  quantitative  analysis  of  the  anomaly  attenuation  is  illus¬ 
trated  in  Figure  5.  This  salt  dome  is  comparable  to  the  Big  Hill  and  High 

SLIDE  #5:  Simulated  analysis 


Island  domes  and  is  generally  typical  of  salt  domes  with  shallow  caprock  in 
the  Gulf  Coast,  U.S.A.  The  density  contrasts  were  taken  from  Nettleton  (1976 
Figure  8-14),  The  calculated  gravity  anomalies,  with  and  without  caprock,  at 
ground  level  and  at  "flight  elevation"  of  1000  feet  (305  m),  are  plotted  in 
solid  lines.  The  attenuation  of  the  caprock  anomaly  is  moderate.  However, 
for  comparison  with  HGMS  gravity  measurements  the  "upward  continued"  ground 
data  must  also  be  filtered.  The  filter  used  by  Carson  had  a  time-width  of 
60  seconds  and  its  width  is  directly  proportional  to  flight  speed.  If  the 
filter  width  is  equal  to  or  larger  than  the  width  of  the  anomaly  the  smooth¬ 
ing  will  be  great.  On  the  other  hand,  if  the  filter  width  is  less  than  the 
width  of  the  anomaly  the  smoothing  will  be  small.  Thus,  we  see  on  slide  5 
that  an  assumed  flight  speed  of  50  knots*  greatly  attenuates  the  sharp  cap¬ 
rock  anomaly  and  spreads  it  over  almost  the  entire  anomaly  (xxxxx).  For  an 
assumed  flight  speed  of  25  knots,  however,  the  filter  width  is  only  half  as 
wide  and  the  attenuation  effect  for  a  caprock  anomaly  of  this  magnitude  is 
very  much  less  ( . ) . 

If  high  frequency,  sharp  anomalies  are  of  interest,  slow  flight  speed  is 
indicated.  The  flight  speed  for  the  1981  demonstration  survey  was  too  fast 
to  record  the  complete  caprock  anomalies  but  the  evidence  of  shallow  positive 
influence  within  the  well-defined  salt  dome  minima  in  Figure  3  are  clearly 
recognizable  to  an  experienced  observer. 

Gravity  gradient  measurements,  as  mentioned  at  the  beginning,  are  beyond 
the  scope  of  this  study.  This  completes  the  overview  portion  of  our  paper. 

*Filter  width  approximately  equal  to  the  width  of  the  caprock  anomaly. 
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COMPARISON  OF  GRAVITY  AND  GRADIENT  ANOMALIES 

We  turn  now  to  the  comparative  study  of  gravity  and  gradient  anomalies. 
Gravity  and  gradient  anomalies  are  quite  different  in  both  magnitude  and  be¬ 
havior  in  space.  Figure  6  compares  gravity  and  several  gradient  components 

SLIDE  Gravity  gradient  comparison 

from  a  horizontal  cylinder  of  infinite  length  at  increasing  altitudes.  In 
the  present  context,  the  point  of  interest  on  this  figure  is  that  the  gradi¬ 
ents,  which  are  seen  to  have  comparable  magnitude  at  ground  level,  attenuate 
with  altitude  an  order  of  magnitude  faster  than  the  gravity  anomalies.  This 
figure  is  from  an  early  paper  (Hammer,  1971)  which  analyzed  the  behavior  of 
gradients  from  a  theoretical  point  of  view. 

Figure  7  is  plotted  on  a  semi-log  scale  and  compares  the  rate  of  attenua¬ 
tion  with  altitude  of  the  gravity  and  vertical  gradient  anomalies  for  a  spheri¬ 
cal  mass.  Gravitv  is  shown  by  solid  curves  and  right  ordinate  scale;  vertical 
SLIDE  :  Gravity  and  vertical  gradient  for  a  spherical  mass 
gradient  by  dashed  curves  and  left  ordinate.  Three  sets  of  curves  represent 
different  depths  (D)  of  the  anomalous  mass.  For  all  cases  the  gravity  anom¬ 
aly  at  ground  level  Lg(D)  is  10  mga 1 . 

If  we  assign  resolution  limits  to  the  gradient  and  gravity  measurements, 
these  curves  define  the  maximum  altitude  at  which  the  anomalies  can  be  detect¬ 
ed.  For  example,  on  the  first  set  of  curves  (with  100  m  depth  to  the  center 
of  mass)  the  magnitude  of  the  gradient  anomaly  drops  below  one  (1)  EBtvBs 
unit  above  altitude  1200  m;  the  gravity  anomaly  drops  below  one  (1)  milligal 
above  altitude  300  m.  The  other  two  sets  of  curves  show  data  for  greater 
depths  and  indicate  higher  limiting  altitudes. 


As  mentioned,  this  represents  a  gravity  anomaly  at  ground  level  of 
10  mgal.  For  larger  (or  smaller)  gravity  anomalies  the  ordinate  magnitudes 
are  to  be  multiplied  by  the  ratio. 

Figure  8  shows  a  similar  analysis  for  a  two-dimensional  gravity  anomaly 
SLIDE  #8:  Gravity  and  vertical  gradient  for  a  horizontal  cylinder 
for  a  horizontal  cylinder  of  infinite  length.  The  limiting  altitudes  of 
detection  for  this  case  are  considerably  larger.  The  maximum  altitudes  at 
which  the  gradient  and  gravity  anomalies  in  Figures  7  and  8  are  detectable 
(i;g  ^  1  mgal,  V22  ^  iE°)  are  listed  in  Table  1. 

TABLE  1:  Maximum  Altitudes  for  Anomaly  Detection 


I  Sphere  .  Cylinder 

I 


D 

meters 

Gravity 
i  meters 

Gradient  1 
meters 

Gravi t  y 
meters 

Gradient 

meters 

100 

[ 

i  316 

1260 

1000 

3162 

500 

;  1581 

3684  ' 

5000 

70  71 

2000 

1  6325 

9283 

20000 

14145 

This  analysis  suggests  that  lower  flight  altitude  is  advantageous  for 
gravity  measurements.  This  may  be  unfavorable  for  fixed-wing  gravity  opera¬ 
tions.  Jordan  (1978)  has  shown  that  low  flight  altitudes  have  other  advan¬ 
tages  as  well. 

This  completes  the  technical  portion  of  our  paper. 


PRESENT  STATUS 


Finally,  it  may  be  of  interest  to  point  out  that  very  extensive  coverage 
of  HGMS  surveys  has  been  accomplished  in  the  past  four  years  (Hammer,  1983). 
Nearly  200,000  km  of  gravity  line  data  have  been  observed  in  20  survey  areas, 

16  of  which  were  overseas  and  4  were  in  the  United  States  and  Alaska.  The 
normal  production  speed  is  50  knots,  but  any  speed  up  to  100  knots  (185  km/hr) 
is  possible.  The  probable  error  of  the  data  is  of  the  order  of  1  mgal  and  the 
anomaly  resolution  limit  is  of  the  order  of  3  km.  The  helicopter  is  capable  of 
flights  at  altitudes  up  to  14,000  feet  (4300  m)  above  any  type  of  terrain  and 
environmentally  sensitive  access.  Duration  of  flights  is  up  to  7  hours.  Nor¬ 
mal  output  per  month  is  about  3,000  miles  (5000  km)  of  line  data. 

A  very  large  project,  which  has  been  in  progress  for  some  time  now,  is 
outlined  on  Figure  9.  If  you  would  like  to  escape  this  cold  winter  weather 

SLIDE  if 9:  Bahamas  project 

in  the  sunny  Bahamas  you  might  like  to  get  involved  with  this  project. 
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CONCLUSION 

In  conclusion:  the  basic  purpose  of  this  paper  is  to  emphasize  to  this 
group  that  the  anomaly  attenuation  of  gravity  gradient  components  is  an  order 
of  magnitude  greater  than  of  the  corresponding  gravitv  anomalies. 

Finally,  I  would  like  to  pose  a  fundamental  dilemma.  The  question  is: 

To  achieve  the  basic  objectives  of  this  symposium  what  is  the  best  way  to 
proceed’ 

(a)  Should  we  continue  the  long-range  development  of  the  theory  and  instru¬ 
mentation  of  the  gradient  measurements;  or 

(b)  Undertake  a  world-wide  gravity  mapping  program  with  the  precise  and 
rapid  HGMS  which  is  now  available? 


**-***pT  5***** 
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FIGURE  CAPTIONS 


Figure  1 

Photograph  of  Carson  HGMS  in  flight.  Note  trailing  magnetometer  at 
I  over  right . 

Figure  2 

Flight  pattern  of  HGMS  survey.  Line  spacing  1  km,  flight  altitude 
1000  ft  (305  m),  and  flight  speed  50  knots. 

Figure  3 

Demonstration  gravity  survey  in  1981  over  known  salt  domes  in  Texas  Gulf 
Coast.  Flight  elevation  1000  ft  (305  m),  flight  speed  45  knots,  and  contour 
i nt erva  1  1  '2  mga 1 . 

Figure  4 

Profile  comparison  of  HGMS  and  conventional  ground  gravity  data  from 
demonstration  survey  in  Figure  3. 

Figure  5 

Analysis  and  comparison  of  simulated  HGMS  and  conventional  ground  grav- 
>.ty  data  relative  to  data  in  Figure  4. 

Figure  6 

Theoretical  comparison  of  magnitude  and  attenuation  of  all  gradient  com¬ 
ponents  and  gravity  anomaly  for  a  horizontal  cylinder  of  infinite  length. 
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Figure  7 

Comparative  attenuation  of  gravity  and  vertical  gradient  anomalies  over 
a  spherical  mass.  Note  semilog  ordinate  which  indicates  that  gradients  at¬ 
tenuate  an  order  of  magnitude  faster  than  gravity.  See  also  Table  1  in  the 
text  . 

Figure  8 

Same  as  Figure  7  but  for  a  two-dimensional  mass  -  a  horizontal  cylinder 
of  infinite  length. 

Figure  9 

Layout  of  HGMS  gravity  survey  now  in  progress  in  the  Bahamas. 

Figure  10  (extra) 

■Attenuation  of  airborne  terrain  correction  simulated  by  an  assembly  of 
two-dimensional  mass  elements. 
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INTEGRAL  FORMULAS  RELATING  GRAVITY 
GRADIENTS  TO  VERTICAL  DEFLECTIONS 


ALAN  H.  ZORN 


Dynamics  Research  Corporation 
60  Concord  Street 
Wilmington,  Massachusetts  01887 


A  number  of  integral  formulas  relating  gravity  gradients  to  vertical 
deflections  are  presented.  The  formulas  are  derived  from  standard  geodesy 
theory  and  are  similar  to  the  Stokes*  and  Vening-Meinesz  integral  formulas. 
The  deterministic  formulas  provide  a  physical  explanation  of  certain 
properties  common  to  various  self-consistent  stochastic  gravity  models. 

For  example,  the  high  degree  of  correlation  between  the  gravity  gradient, 

T  ,  and  the  vertical  deflection,  T  ,  in  most  gravity  models  can  be 

XZ  X 

physically  justified  based  on  a  deterministic  integral  formula.  The 
integral  formulas  may  lead  to  insight  into  the  construction  of  gradiometer 
filters  to  estimate  vertical  deflections  which  are  insensitive  to  stochastic 
gravity  models.  Implementation  considerations  involving  discretization 
and  truncation  of  the  integral  formulas  are  explored  using  surface  gravity 
data  over  the  Blake  Escarpment. 
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TRANSFORMATIONS  OPERATE  IN  THE  LEVEL  x-y  PLANE  ONLY 
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(Ci:!7I  IVIVHON)  NOIXOM'I.-JMd  IVJUHMA 


IS  RADIUS  OF  CIRCLE  BOUNDING  REGION  OF  INTEGRATION 
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SYNOPSIS 


REAL-TIME  VERTICAL  DEFLECTION  ESTIMATION 
USING  GRADIENT  DATA 


WALTER  K.  FELDMAN 
Sperry  Corporation/ Electronic  Systems 
Great  Neck.  New  York  11020 
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Techniques/algorithms  for  making  real-time  estimates  of  vertical 
deflections  from  gravity  gradient  measurements  are  discussed.  First, 
techniques/algorithms  used  to  condition  sensor  data  for  effects  of 
high-frequency  noise,  pressure  sensitivity,  self-gradients,  bias/trend, 
and  carouselling  are  outlined.  Plots  of  typical  gradient  data  before  and 
after  preprocessing  are  shown  to  illustrate  the  effectiveness  of  the 
preprocessing.  Second,  a  batch  filter  formulation  for  estimating 
deflections  from  gradient  and  SEASAT  map  data  is  illustrated.  Finally, 
filter  results  are  demonstrated  by  plots  of  filter  estimates  of  vertical 
deflection  vs.  reference  values  for  several  ship  tracks. 
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OVERVIEW:  PART  I 
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TABLE  1.  NED  ERRORS  DUE  TO  INSTRUMENT  ERRORS 
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PROVIDES  ANOMALOUS  GRADIENTS  DIRECTLY 

PROVIDES  ESTIMATION  OF  RESIDUAL  NOISE  FOR  MONITORING 
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FIGURE  2. 


CURE  2A  NORTH- 
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FIGURE  2B  NORTH-DOWN  AND  EAST-DOWN 
aDIENTS  BEFORE  AND  AFTER  PREPROCESSING 
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ABSTRACT 

One  of  the  purposes  of  NASA/Goddard  Space  Flight  Center's  proposed 
Geopotential  Research  Mission  is  to  map  the  gravity  field  of  the  Earth  to  1-2 
mgal  in  41,253  1  degree  x  1  degree  blocks  using  range  rate  tracking  data 
between  two  polar  low  orbiting  satellites  in  the  same  orbit.  Because  there  are 
41,253  blocks  in  the  gravity  map,  it  takes,  in  general,  at  least  9.4x10^^ 
real  floating  point  operations  (flop)  to  solve  the  least  square  fit  (LSQF) 
matrix  equation  used  to  determine  the  gravity  potential  from  the  tracking 
data.  In  this  paper,  by  studying  the  synmetry  of  an  operator  which  generates 
the  range  rate  acceleration  between  the  satellites  from  the  along  track 
gravity  field,  it  is  shown  that  the  LSQF  cross  correlation  matrix  (T-matrix) 
formed  from  a  complex  spherical  harmonic  {Y|j^)  expansion  of  the  gravity 
notential  is,  in  the  circular  orbit  infinite  trajectory  length  approximation, 
diagonal  in  n.  Further,  for  the  actual  mission,  an  iterative  method  which 
utilizes  the  inverse  of  the  m-diagonal  portion  of  the  T-matrix,  is  shown  to 
produce  the  32758  coefficients  of  the  Y^^  model  of  the  gravity  potential  to  64 
bit  precision  in  less  than  IxlO flop.  A  second  iterative  method  is  shown  to 
oroduce  an  approximate  covariance  matrix  for  the  coefficients  in  less  than 
another  5x10“  flop.  This  paper  also  demonstrates  a  method  which  eliminates 
the  problem  of  aliasing  when  attempting  to  solve  the  full  problem 
sequentially,  first,  by  solving  for  the  approximate  trajectories  of  the 
satellites  with  a  reduced  model  of  the  gravity  fit  Id,  and  then,  by  solving  for 
a  full  model  of  the  gravity  field  using  the  approximate  trajectories. 


1.  INTRODUCTION 


One  of  the  purposes  of  Goddard  Space  Flight  Center's  (GSFC)  proposed 
Geopotential  Research  Mission  (GRM)  (Smith,  et.  al.,  1982)  is  to  map  the 
gravity  field  of  the  Earth  to  1-2  mgal  in  41,253  1  degree  x  1  degree  blocks 
using  range  rate  tracking  data  between  two  polar  low  orbiting  satellites  in 
the  same  orbit  (low-low  range  rate  tracking).  Using  worst  case  analysis,  both 
GSFC  and  the  Johns  Hopkins  University  Applied  Physics  Laboratory  (APL),  which 
has  performed  much  of  the  design  and  study  work  for  the  proposed  mission,  have 
concluded  that  the  data  analysis  required  to  turn  the  range  rate  data  into  the 
gravity  map  will  consume  very  large  amounts  of  computer  time  even  on  the 
fastest  computers  that  exist  today  (Wei f fenbach,  1983;  vonBun,  1983).  One  of 
the  ^nst  time  consuming  parts  of  the  data  analysis  is  in  solving  the  least 
square  'it  (LSQF)  matrix  equation  used  to  determine  the  gravity  potential  from 
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the  tracking  data  because  of  the  4K253  blocks  'n  the  gravity  map.  It  is 
later  shown  that  the  solution  of  this  equation  using  Gaussian  elimination 
techniques  (See  Section  5  and  Appendix  B. )  will  take,  in  general,  at  least 
9.4x10  ^^real  floating  point  operations  (flop). 

By  examining  the  symmetry  of  an  operator  which  generates  the  range 
acceleration  from  the  along  track  gravity  field,  this  paper  demonstrates  that 
the  LSQF  cross  correlation  matrix  (T-matrix)  formed  from  a  complex  spherical 
harmonic  expansion  of  the  gravity  potential  is,  in  the  circular  orbit 

infinite  trajectory  length  approximation,  diagonal  in  m  Further,  for  the 
actual  mission,  an  iterative  method  which  utilizes  the  inverse  of  the 
m-di agonal  portion  of  the  T-matrix,  is  shown  to  produce  the  32758  coefficients 
of  the  Yj  model  of  the  gravity  potential  in  less  than  1x10^^  flop.  A  second 
i terative  method  is  shown  to  produce  an  approximate  covariance  matrix  for  the 
coefficients  in  less  than  another  5x10^^  flop. 

This  paper  also  demonstrates  a  method  which  eliminates  the  problem  of 
aliasing  when  attempting  to  solve  the  full  problem  sequentially  by  solving  for 
the  approximate  trajectories  of  the  satellites,  and  then,  first,  by  solving 
for  the  full  gravity  potential  using  the  approximate  trajectories. 


2.  GRM  GRAVITY  FIELD  NORMAL  EQUATIONS 


As  shown  in  Figure  1,  the  2  GRM  satellites  will  be  launched  in  the 
same  160  Km  polar  nearly  circular  orbit  with  an  along  track  separation  of 
about  300  Km  (Smith,  et.  al.,  1982).  On  board  navigation  and  control  systems 
effectively  make  the  motion  of  both  satellites  drag  free  so  changes  in 
satellite  motion  for  times  longer  than  a  few  seconds  are  solely  due  to  the 
gravity  field  of  the  Earth  (Ibid.).  Thus,  the  motion  of  each  satellite  can  be 
described  by: 


f  =  -?V(r  ,t)  (2.1) 

1  1 

where  r^,  satellite  i's  center  of  mass  position,  and  v,  the  gradient,  are  in 
units  of  R the  average  orbital  radius.  The  coordinate  system  for  this 
equation  is  an  approximately  inertial  geocentric  system,  so  the  gravity 
potential  of  the  Earth, V,  is  a  function  of  time  due  to  the  Earth's  rotation. 
From  this  point  on,  the  time  dependence  of  V  will  not  be  explicitly  written. 

To  determine  the  Earth's  gravity  field,  GRM  will  use  the  observable  p, 
the  magnitude  of  the  range  vector  between  the  center  of  mass  of  the  satellites 
gi ven  by; 


P  =  (2.2) 

Since  the  satellites  are  in  the  same  nearly  circular  orbit  and  p  is  small 
compared  with  the  circumference  of  the  orbit,  p  is  essentially  the  along  track 
component  of  ~  we  can  write: 

p(t)  =  -  VjVir^)  +  (2.3) 

where  Ts  is  the  component  of  the  gradient  along  ds,  the  differential  along 
track  vector  from  satellite  1  to  2  in  units  of  R h.  We  can  rewrite  this 
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in  a  more  compact  form  using  a  position  difference  operator,  P(^),  defined  by: 

P{^)f(‘j^)  =  f(r+p)  -  f(r).  (2.4) 

Using  this  operator,  (2.3)  becomes  a  normal  equation  for  p  in  terms  of  V  and 
the  operator  P! 


p(t)  =  -  PV^ir)  12.5) 

where  P  will  be  implicitly  P{^)  unless  the  functional  dependence  is  otherwise 
specified.  A  normal  equation  for  p,  the  range  rate  observable,  can  be 
generated  from  (2.5)  by  integrating  the  equation. 


3.  CIRCULAR  ORBIT  APPROXIMATION 

If  we  assume  that  the  orbit  is  circular,  (2.5)  can  be  written  in  a 
form  which  will  prove  extremely  useful  later. 

In  a  polar  coordinate  system,  (r,  e,  <|)),  whose  north  polar  axis 
coincides  with^the  North  Pole  of  the  Earth,  the  along  track  diferential 
coordinate,  ds,  is  parallel  to  the  0  coordinate.  If  the  descending  node  of  the 
orbit,  (fn  is  0  deg  -  180  deg,  ds  =  d0.  If  the  descending  node  of  the  orbit  is 
180  deg  -  360  deg,  ds  =  -de.  Thus,  since  r  =1  at  the  orbit  height: 

V  *  ±  i  (3.1) 

5  30 

where  the  plus  or  minus  depends  on  the  value  of  as  given  above. 

Using  (3.1)  and  integrating  (2.5)  with  respect  to  t,  we  obtain  an 
exDlicit  form  for  the  normal  equation  for  p  in  the  circular  orbit 
approximation  given  by: 


p(t)  =  -  PV/v  (3.2) 

where  v  is  the  satellite  orbital  speed.  To  obtain  (3.2),  we  have  used  the  fact 
that: 


dt  =  ds/v  =  *d0/v  (3.3) 

and  that  P  commutes  with  the  integration. 

We  can  also  derive  an  explicit  operator  expression  for  P  by 
considering  the  properties  of  the  displacement  operator,  D,  defined  by: 

D(p)f(f)  =  f(f+p).  (3.4) 

To  derive  the  properties  of  D,  first  consider  the  operator  for  a  differential 
displacement.  (This  is  a  standard  derivation  found  in  many  texts,  i.e., 
Tinkham,  1964;  Messiah,  1966.)  For  a  differential  displacement  in  the  along 
track  direction,  dl,  D  is  given  by: 

0(ds)f(f)  =  ds-v  f(r)  +  f(r)^ 


(3.5) 
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so; 


D(dl)  =  1  +  dsV^.  (3.6) 

If  we  can  assume  all  the  d?  are  along  the  direction  of  p  (This  true  for  p«l, 
the  case  here.),  then  we  can  write  p  =  nds.  Thus  D(p)  can  be  written  as 
(D(ds))  yielding: 

pv 

0(p)  =  Tim  (1  +  pV^/n)'^  =  e  ^  (3.7) 

n-x« 


Using  (3.1),  we  obtain; 


0(p) 


(3.8) 


(If  p  were  not  small,  in  the  above  expression,  p  would  have  to  be  replaced  by 
the  arc  length  between  the  end  points  of  p. ) 


But  P  can  be  written  as: 


P(p)  =  D(?)  -  1. 

(3.9) 

Therefore: 

+  D— 

P(?)  =  e'  -1. 

(3.10) 

4.  SEQUENTIAL  LEAST  SQUARES  FITTING 


The  standard  technique  of  deriving  V  from  (2.5)  or  (3.2)  using  the 
range  rate  data  is  to  represent  V  by  a  set  of  known  basis  functions  times 
unknown  coefficients  and  to  use  least  squares  fitting  techniques  to  determine 
the  coefficients  and  the  trajectory,  r(t).  Since  our  primary  concern  is  V  and 
not  r(t),  one  can  attempt  to  solve  the  full  problem  sequentially,  first,  by 
solving  for  the  approximate  trajectories  of  the  satellites  with  a  reduced 
model  of  the  gravity  field,  and  second,  by  solving  for  a  full  model  of  the 
gravity  field  using  the  approximate  trajectories.  Let  us  consider  the 
consequences  of  using  this  approach. 

For  this  discussion,  let  our  normal  equation  be  represented  as: 

p(t)  =za.Pf.(f)  +  z:(b  .+c  .)Pg  .(r)  (4.1) 

where  we  want  to  determine  the  cofficients  a^  and  b^+c^  given  range  rate  data 
taken  over  some  unknown  trajectory  ^(t)  and  the  known  sets  of  functions  f^  and 
g^.  Assume  also  that,  through  least  squares  fitting  over  the  same  trajectory, 
possibly  using  ancilliary  data,  we  have  found  a  solution; 

Po(t)  =  Eb  .Pg  .(^  )  (4.2) 

which  has  yielded  the  values  b^  and  the  approximate  orbit  vector  r ^[t) .  (The 
fit  may  involve  each  satellite  individually  and  nay  use  a  priori  b^'s.) 
Defining  ip  =  -  Pq  and  subtracting  (4.2)  from  (4.1),  yields: 


5 


iipit)  =  Ea.Pf.(r)  +  z(b>c.)Pg.(r)  -  zb^.Pg  .(r^).  (4.3) 

i  i 

Further  assume^that  the  functions  f  .(r)  and  g^-(r)  can  be  approximated  by 
f;{fQ)  and  gi(rQ).  (f  and  g  don't  change  appreciably  over  the  distance 

Ir-r  1.)  This  turns  (4.3)  into: 

0 

Ap(t)  =  za^Pf.  (r^)  +  rc^-Pg^-C^’p)  (4.4) 

It  is  necessary  to  include  the  set  of  c-coefficients  in  the  problem 
because  the  coefficients  of  the  g-fuctions  derived  from  a  least  square  fit  of 
(4.1)  will  not,  in  general,  be  the  same  as  the  the  coefficients  derived  from  a 
least  squares  fit  which  doesn't  include  the  f-functions;  to  obtain  the  least 
squares  condition  without  the  f-functions,  the  g-f unction  must  alias  the 
effects  described  by  the  f-function  through  alteration  of  the  b-cofficients 
from  the  "true"  value  obtained  if  the  full  least  squares  fit  were  performed. 
(This  is  also  true  if  a  priori  b-coefficients  used  in  the  1st  fit,  were 
determined  from  past  fits.)  Thus  by  including  the  the  c-  coeff»cien+s  H\e 
g-functions  in  a  second  least  squares  fit  on  (4.4),  we  can  correct  our 
original  b-coefficients  to  obtain  better  values  and  we  can  ensure  that  the 
previous  fit  does  not  cause  aliasing  problems  in  the  a-coefficients.  Since  the 
coefficients  of  the  g-functions^are  being  refit,  the  only  importance  of  the 
first  fit  is  to  produce  a  good  r^it). 

If  we  were  interested  in  a  further  refinement  in  the  trajectory,  we 
could  apply  the  approach  used  to  obtain  (4.4)  to  generate  an  expression  for  a 

correction  to  the  trajectory  given  by: 

«  « 

r(t)  =  -  ra^V(vfj(rj)-Z(b.+c  .)v(vg.(»*g))  (4.5) 

This  expression  could  be  integrated  twice  to  generate  a  better  estimate  of 
r(t)  after  the  a^'s  ,  b^'s,  and  c.'s  were  determined.  This  better  r(t)  could 
then  be  used  to  obtain  an  improvea  data  set  for  generating  AP(t)  which  in  turn 
could  be  used  to  obtain  improved  a  -  and  c^-  values,  etc.  Thus  (4.4)  and  (4.5) 
can  be  used  together  in  an  iterative  chain  to  refine  the  results  obtained 
before  iteration. 

4.2  USING  LEAST  SQUARE  FITTING  TO  DETERMINE  THE  GRAVITY  FIELD 

If  we  now  assume  the  g-functions  are  to  be  part  of  the  f-functions  and 
set  Ap  =  y,  we  can  write  (4.4)  in  the  general  form: 

y(t)  =  l  a  ^Pf  ^  (rQ)  (4.6) 

Let  us  allow  the  fit  functions  and  the  coefficients  to  be  complex  even  though 
y  is  real.  (This  will  lead  to  simpler  notation  later.)  The  sum  of  the  squares 
or  is  then  given  by: 


where  w’(ro)  is  a  weighting  factor  t^  be  determined  later  and  where  the 
integration  is  over  the  trajectory,  r^. 


(4.7) 
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If  the  elapsed  time  for  the  trajectory  is  long  compared  with  24 
hours,  the  rotation  period  of  the  Earth,  the  trajectory  will  pass  over  or 
nearly  over  every  point  on  Earth  with  a  probability  which  is  only  function  of 
0  (infinite  trajectory  limit).  (The  orbit  period  of  the  satellite  must  also 
not  be  an  exact  harmonic  of  the  Earth's  rotation  period,  the  case  here.)  In 
the  infinite  trajectory  limit,  the  orbit  also  passes  over  each  point  an  equal 
number  of  times  with  ds=+de  and  ds=-d0.  In  this  limit  and  assuming  the 
circular  orbit  approximation,  (4.7)  thus  becomes: 


=  E  /  dfi  w(n)(y-Ea*Pf*(fi)){y-Za.Pf.(fi))  (4.3) 

±  i  ''  i 

where  we  have  replaced  the  integral  over  the  trajectory  with  a  sum  over  the 
two  possible  directions  of  ds  in  P  of  an  integral  over  a  geostationary 
spherical  surface  of  radius  times  an  orbit  densi'^  function  absorbed  into 
w(/i).  In  (4.8),  notice  that  y  is  now  a  function  of  r,  not  t,  and  the  implicit 
t  dependence  in  Y  has  disappeared,  y  must  now  be  interepreted  as  the  average 
value  of  ^  over  the  earth  fixed  point  {6,41). The  new  weighting  function  w  (n)  is 
given  in  terms  of  the  old  one  by: 

w(  i^)=n(  0)w'  (r)  (4.9) 

0 

where  n(6)  is  the  density  of  trajectory  points  over  the  Earth  region  (6,41)  given 
by: 

n(  e)  =  Nrev/{ffsin0)  (4.10) 

where  Nrev  is  the  number  of  revolutions  in  the  trajectory.  (4.3)  is  not  a 
practical  for  performing  the  actual  least  square  fit,  but  will  be  useful 
for  the  theoretical  analysis  of  the  least  square  fit  problem  to  follow. 

From  (4.3),  least  square  fit  equations  can  be  developed  in  the  usual 
way  (Wolberg,  1967)  by  taking  the  derivatives  of  x^with  respect  to  the 
a-coefficients.  This  results  in  the  following  matrix  equation: 

L  =  TA  (4.11) 


where: 


A,  =  a. 


Ly  =  ( l/2)E/dfiwyPf.(e,4) 
T,-/  =  (l/2)E;dfiwP*f*  Pf. 


(4.12) 

(4.13) 


(4.14) 


L  and  A  are  column  vectors  and  T  is  a  hermetian  matrix  whose  inverse  provides 
the  solution: 


A  =  T"^  L 


(4.15) 


given  by: 


The  covariance  matrix  C-jj  =  {<..>  is  an  ensemble  average)  is 


C  =  t'^xVf 


r  u  -  I  A  /r 

e  O  /Aj-  ,  s  ^  Ae  e  F  A  •  fk,  . 

'  *•  w  I  If  e  • 


(4.16) 
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where  F  is  the  number  of  degrees  of  freedom,  the  number  of  independent 
measurements  minus  the  number  of  a-coefficients.  We  can  determine  by 
expanding  (4.8)  and  using  (4.15)  to  yield: 

=  y"-A^L  =  72  -  L+A  (4.17) 

+  * 

where  A  L  =  EA^L^-,  etc,  and  where: 

■  iE/w(e)  y^dft  (4.18) 

+ 

If  the  0^  (4.7)  is  used  instead  of  that  of  (4.8),  the  equations 
remain  the  same  except  that  the  sum  plus  integral  over  the  spherical  surface 
goes  back  to  an  integral  over  the  trajectory  and  w  is  replaced  by  w' 


5.  SOLVING  THE  MATRIX  EQUATIONS 


THE  SCALE  OF  THE  OVERALL  PROBLEM 


The  basic  GRM  gravity  data  analysis  problem  is  to  find  the  orbit 
approximation,  r'^(t),  to  solve  the  matrix  equation  (4.11)  in  terms  of 
functions  at  the  orbit  height,  and  to  analytically  extend  the  solution  at 
orbit  height  to  the  geoid.  The  analytic  extension  problem  can  be  solved 
without  large  computations  by  using  basis  functions  for  V  which  are  solutions 
to  Laplace's  equation.  We  will  not  concern  ourselves  with  the  orbit 
computation  problem  here.  We  will  concern  ourselves  here  with  the  remaining 
problem,  the  solution  of  the  matrix  equation. 


There  are  two  aspects  of  the  matrix  problem,  forming  the  matrix 
elements  and  inverting  T.  The  size  of  these  tasks  are  determined  by  the  number 
of  basis  functions,  N,  and  the  number  of  range  rate  data  points  collected,  n. 

N  is  determined  by  the  number  of  blocks  in  the  gravity  map.  In  a  coordinate 
lattice  which  minimizes  the  number  of  blocks  (regular  square  spherical  shell 
lattice),  the  1  degree  resolution  requirement  of  GRM  (Smith,  et.  al.,  1983) 
translates  into  41,253  such  blocks.  If  localized  functions  are  used  to  make 
the  map,  one  function  to  a  map  block,  then  there  will  be  41,253  basis 
functions.  If  complex  spherical  harmonics,Yj^jj^,  are  used  to  the  point  where  one 
half  wavelength  of  the  function  ( n/ z)  equals  the  block  size,  there  will  be 
spherical  harmonics  to  orderil„=180  and  there  will  be  ( ((y-l)^-3)=32758  real 
basis  functions.  For  a  6  montn  mission  with  1  data  point  per  4  seconds 
( Ibid. ),  n  is  3.9x10®  . 

Inverting  an  NxN  matrix  using  Gaussian  elimination  techniques  requires 
approximately  (4/3)N^  real  floating  point  operations  (flop)  (See  Appendix  8.) 
and  requires  access  to  the  matrix  elements  N  times.  Use  of  caching  techniques 
as  planned  by  GSFC  in  the  proposed  mission  (Wei ffenbach,  1983)  and  as 
described  by  others  (Duff,  1981)  can  reduce  the  number  of  of  disk  accesses  to 
the  point  where  a  large  processor  can  run  at  10  to  20  percent  of  full 
processing  speed  (Briggs,  1984).  This  means  that  with  a  Cyborg  type  computer, 
one  can  plan  on  having  10-20  Mflop/s  of  effective  speed  (Ibid.).  Using  this 
number  for  the  processing  speed  and  41,253  for  N,  we  obtain  that  9.4x10^^ 
flops  and  1300-2600  hrs  of  computing  time  required  for  inversion  of  T. 


Forming  the  matrix  elements  for  the  matrix  equation  is  dominated  by 
the  calculation  of  T.  Naively  approaching  the  problem,  it  would  seem  that 
there  are  on  the  order  of  nxN  =  7x10  flop  required  to  calculate  the 
elements  of  T,  based  on  calculating  H^/2  basis  functions  for  n  data  points. 

Even  at  80  Mflop,  this  would  require  23,000  hrs  to  compute  the  elements. 
However,  since  the  basis  functions  only  change  for  lengths  on  the  order  of  a 
block  size,  one  does  not  have  to  calculate  the  basis  functions  for  each  point. 
APL  has  developed  a  scheme  to  do  this  in  approximately  5000  hrs  of  computer 
time  (Wei ffenbach,  1980).  The  problem  of  calculating  the  natrix  elements, 
though  clearly  the  dominant  problem,  will  not  be  addressed  here. 

One  approach  to  finding  a  way  of  reducing  the  computation  time 
associated  with  inverting  T  is  to  investigate  the  syime tries  and  near 
symmetries  associated  with  the  P  operator  and  various  sets  of  basis  functions 
to  see  if  the  matrix  can  be  broken  down  or  approximately  broken  down  into 
several  irreduceable  submatrices  (diagonal  blocks  of  submatrices)  of  lower 
dimensionality  (Tinkham,  1964).  These  can  be  inverted  with  less  computation 
time  because  the  number  of  operations  goes  as  such  a  high  power  (N*)  of  the 
dimensionali ty. 

If  a  symmetry  is  only  approximate,  it  will  create  small  terms  instead 
of  zeros  off  the  irreduceable  submatrices.  If  a  residual  matrix  equation  is 
created  by  operating  on  the  original  matrix  equation,  TA  =  L,  with  the  inverse 
of  the  approximately  i rreproduceable  submatrices  of  T,  it  is  shown  in  Appendix 
B  that  the  residual  equation  can  be  solved  iteratively  if  the  residual 
T-matrix  is  diagonally  dominant.  If  this  is  the  case.  Appendix  B  also  shows 
that  the  residual  equation  can  be  solved  for  A  iteratively  in  less  than  1x10^^ 
flop  (Appendix  B.4)  and  that  an  approximate  T'^  for  determining  the  covariance 
can  be  generated  in  another  5x10^^  flop  (Appendix  B.6). 

Investigations  by  the  author  using  the  form  of  the  P  operator  derived 
in  Section  3  were  carried  out  using  point  anomaly  basis  functions  and 
spherical  harmonic  basis  functions.  The  point  anomaly  functions  were  generated 
from  the  gravity  field  of  point  masses  at  the  center  of  square  spherical  shell 
lattice  blocks  at  a  depth  below  the  surface  equal  to  the  width  of  the  block. 

The  T-matrix  generated  from  these  anomaly  functions  was  highly  diagonal  but 
unfortunately  was  not  diagonally  dominant.  Also  unfortunate  was  the  lack  of 
correlation  in  the  large  off  diagonal  terms  (nearest  neighbors)  which  may  cause 
back-filling  problems  (See  Appendix  8.1.).  (This  is  an  essential  feature  of 
trying  to  index  a  two-dimensional  surface  with  a  one  dimensional  matrix  index 
-  most  nearest  neighbors  are  not  nearest  in  index  number.)  To  show  that  the 
matrix  does  or  does  not  have  back-filling  problems  will  require  further 
invest!  gation. 

The  investigation  using  spherical  harmonic  basis  functions  proved  so 
successful  in  demonstrating  that  there  is  an  approximate  symmetry  which  will 
greatly  reduce  the  computation  time  required  for  inverting  T  that  only  this 
approach  will  be  presented  here.  The  next  section  presents  the  approach  based 
on  these  functions. 


f 
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6.  CONSEQUENCES  OF  USING  SPHERICAL  HARMONIC  BASIS  FUNCTIONS 


Let  us  consider  the  consequences  of  using  spherical  harmonics  as  the 
basis  functions  for  V.  To  simplify  writing  the  equations  and  for  other  reasons 
which  shall  become  apparent  later  in  this  section,  let  us  use  the  complex 
spherical  harmonics,  Y  (e,^),  with  complex  coefficients  a^^  rather  than  the 
real  spherical  harmonies,  and  used  normally  in  Geodesy 

(Heisman  and  Moritz,  1967).  in  terms  of  a  solution  to  Laplace's  equation 

valid  in  all  regions  outside  a  sphere  of  any  radius  is  given  by  (Jackson, 
1962): 

V(r.e,»)  =  "i  I  Vj„(e.*)  <6-l' 

«,=o  m=-e.  r 

where  we  have  normalized  the  coefficients  so: 


f(1.0,$)  =-  = 


OO  H 

L  E 

Z=o  m=-Z 


Y,  is  given  by  (Messiah,  1966,  Vol.  1,  Appendix  B): 

tc  m  I 


m4) 


(6.2) 


(6.3) 


where  P.^nis  the  associated  Legendre  polynomial  of  order  (^m)  given  by  the 
following  for  positive  m  (Ibid.): 

p,„(x)  =  (6.4) 


£m' 


dx 


For  negative  m,  the  P.  are  defined  by  (Ibid.): 
^  1  m 


Y  =  (-1)  ^  Y* 

£  ,-m  ^  '  SLm 


(6.5) 


As  defined  here,  the  Y^^are  orthonormal  with  respect  to  integrations  on  a 
uni t  sphere  ( Ibi d. ) . 


For  a  spherical  harmonic  basis,  in  the  circular  orbit  approximation, 
the  T-matrix  elements  become: 


T 


t'm'lm 


★  ★ 


P  dn 


(6.6) 


where  we  have  put  a  restriction  on  w(n)  that  it  be  only  a  function  of  9. 
Normally,  the  weighting  function  in  terms  of  the  trajectory  integral,  vv'(rj), 
will  be  a  constant;  since  p  is  approximately  the  same  everywhere  along  the 
trajectory,  we  expect  the  range  rate  data  collected  all  along  the  trajectory 
to  have  about  the  same  noise  content.  This  means  that  normally  w(fi)  is 
proDortional  to  n(9),  so  the  restriction  is  normally  satisfied.  In  the  next 
paragraph,  we  will  show  that  this  restriction,  along  with  the  properties  of  P 

in  the  circular  orbit  approximation,  are  instrumental  in  simplifying  the 
inversion  of  the  T-matrix. 
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The  expression  for  the  P  operator  in  the  circular  orbit  approximation, 
(3.10),  is  written  in  terms  of  the  operator  3/ 3a  Since  P  has  no  <))  dependence, 
it  does  not  effect  the  functional  dependence  of  the  with  respect  to 

(e  Thus,  since  we  assumed  w(j2)  is  only  a  function  of  0,  the  matri){ 

elements  will  be  zero  if  m  is  not  equal  to  m*  because  the  integral  of 
over  360  degrees  is  zero  unless  m=m',  regardless  of  the  0  behavior.  In  other 
words,  the  T-matrix  is  diagonal  in  m.  For  a  given  m,  there  are  only  Yj^m's  with 
with  ii>=|m|,  so  each  is  non-zero  only  for  I  >  !m!  and  )i’>|m|. 

The  fact  that  the  T-matrix  is  diagonal  in  m  greatly  reduces  the  number 
of  conputer  operations  required  to  invert  the  T-matrix.  If  we  subgroup  the 
indices  by  m  first  rather  than  by  l  first,  we  can  form  the  T-matrix  into 
22,  +1  diagonal  m  submatrices.  (2,^  =180  is  the  maximum  value  |m|  can  have.) 
Ea^  m  submatrix  for  m>=l  has  (  2^ -lml+1)  ^  complex  elements  (  (22^  +1) 
elements  for  m=0).  Summing  for  m  =  -2^  to  2^  the  approximately 
(16/3)(  2  -|m(+l)  ^  real  operations  required  to  invert  each  submatrix,  we  find 
that  only  approximately  (8/3)  2^  =2.8x10®  flop  are  needed  to  invert  the  full 
matrix.  Also  of  impact  on  the  computation  time  is  the  fact  that  the  largest 
submatrix  (m=0)  only  requires  260642  real  words  for  storing  the  matrix 
elements  (and  only  64800  words  for  the  next  largest,  etc.),  so  we  can  load 
each  submatrix  to  be  inverted  wholly  into  core  memory  on  a  large  machine.  This 
means  the  machine  can  run  at  virtually  full  speed  to  invert  the  matrix.  So 
assuming  an  80  Mf lop/s  computation  rate  and  a  1  Mbyte/s  disc  transfer  rate, 
one  can  invert  the  matrix  in  66  seconds. 

To  aid  in  further  investigations  of  this  approach,  it  would  be  useful 
to  be  able  to  calculate  the  matrix  elements  of  T  in  the  circular  orbit 
aoproximati on.  Appendix  A  shows  how  to  generate  closed  formulas  for  the  matrix 
elements  without  actually  performing  the  integrations.  It  shows  that  the 
operator  P  in  the  circular  orbit  approximation  can  be  written  in  terms  of  the 
quantum  mechanical  angular  momentum  operators,  L^and  L_  (Messiah,  1966),  which 
have  well  known  properties  when  acting  on  Y  and  yield  linear  combinations  of 
other  Y  ij^'s.  It  also  shows  that  the  orbit  density  function  in  the  circular 
orbit  approximation  times  Y  ,  can  be  expressed  as  linear  combinations  of  other 
Y  j^'s.  Together  these  properties,  along  with  the  orthonormal  property  of  the 
y]|^'s,  allow  T  ,  to  be  written  in  terms  of  sums  over  coefficients  which 
are  functions  of  .’,m,  2' ,m' ,and  p. 


7.  WILL  AN  ITERATIVE  INVERSION  OF  THE  RESIDUAL  T-MATRIX  CONVERGE? 


We  showed  for  (A)  a  circular  orbit  polar  orbit,  (3)  w(r)  =  w(0),  and 
(C)  the  infinite  trajectory  limit  that  the  T-matrix  can  be  inverted  fast 
because  it  is  diagonal  i n  m.  In  Appendix  B,  it  is  shown  that  the  actual  matrix 
equation,  TA  =  L,  can  be  solved  iteratively  for  A  (Appendix  B.4)  and  for  an 
approximate  covariance  (Appendix  B.6)  if  the  residual  T-matrix  formed  by 
operating  on  T  with  the  inverse  of  the  diagonal-m  portion  of  T  is 

diagonally  dominant  (Appendices  3.3  and  B.5).  Let  us  examine  the  violations  of 
the  assumptions  (A),  (B),  and  (C)  that  will  occur  in  the  actual  planned 
mission  to  see  if  the  diagonal  dominance  criterion  is  satisfied  or  not. 
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In  the  infinite  trajectory  limit,  for  the  diagonalization  in  m  to  be 
violated  there  must  be  some  ^dependence  in  P,  w(r^),  or  Vs  (See  equati on2.5. ). 

A  small  eccentricity  in  the  orbit  would  not  cause  problems  because  this  only 
causes  mixing  of  the  radial  and  e  coordinates.  A  deviation  of  the  orbit  from  a 
polar  one  causes  mixing  of  the  <|)  coordinate  into  the  p.?  argument  of  P  and 
into  Vs,  but  in  a  way  which  doesn't  effect  the  diagonalization  asPoilows.  Tl^e 
dependent  part  of  p*v  can  be  written  as,  and  the  <))  dependent  part  of  Vs  is 
proportional  to: 

(p.v)  =  -tp  (0)(i/sin  0  )  (7.1) 

where  p  *(  6)  functionally  describes  the  orbit  with  e  as  the  independent 
variable^  and  where  .the  *  is  .for  the  2  possible  directions  of  ds  at  the  point 
(e,  4i).  Since  ^e  =ime  ,  this  term  does  not  change  the  orthogonality 
of  the  so  the  matrix  elements  off  diagonal  in  m  remain  zero.  Also  a 

deviation  of  the  orbit  from  a  polar  one  does  not  cause  w(5)  to  become  a 
function  of  <p  since  the  density  remains  axially  symmetric  regardless  of  the 
orbit  in  the  infinite  trajectory  limit.  Therefore,  a  deviation  of  the  orbit 
from  a  polar  one  does  not  produce  terms  off  diagonal  in  m  in  the  infinite 
trajectory  limit. 

The  finiteness  of  the  actual  data,  however,  can  introduce  $  dependence 
into  w(i^J.  As  stated  previously,  this  effect  goes  to  zero  in  the  limit  of  an 
infinite  trajectory  length  since  the  phase  relationship  between  the  rotation 
of  the  Earth  and  the  rotation  of  the  spacecraft  around  the  geocentric  fixed 
orbit  go  through  all  possible  relationships  equally.  Since  the  orbital  period 
is  much  less  than  the  Earth's  rotation  period,  the  effect  of  the  finite  data 
length  in  the  limit  of  many  orbit^  is  given,  to  1st  order,  by  the  residual 
assymmetry  in  the  trace  of  the  orbit  as  the  Earth  turns  under  the  orbit.  Since 
the  Earth  turns  under  the  orbit  once  in  24  hours,  the  relative  assymmetry  in 
n(fl)  will  be  given  by  a  maximum  of  (12  hrs)/(6  mo)  =3xl0“\  This  number  is  just 
above  the  level  of  the  order  of  magnitude  criterion  for  diagonal  dominance 
based  on  the  assumption  that  all  the  off  diagonal  elements  in  m  are  equal  as 
given  in  (8.18).  (See  Appendix  B  for  a  more  detailed  discussion.)  We  must, 
thus,  examine  this  effect  in  more  detail. 

The  assymmetry  is  a  periodic  step  type  function  in  with  a  2tt  period 
and  an  amplitude  0.003.  (For  ^a<<t><Tb,  w  is  less  by  0.003  than  w  for  other 
values  of  -)>).  Since  the  <{'  dependence  of  is  given  by  ,  the  relative 

size  of  to  1  will  be  given  by  the  the  absolute  value  of 

(m-m')th  harmonic  amplitude  of  the  fourier  transform  of  a  periodic  step  type 
function.  The  absolute  value  this  harmonic  is  less  than  0.003*2/ (tt !m-m' ! ) 
regardless  of  the  width  of  the  step  (Selby,  1974).  Therefore  the  criterion  for 
convergence  of  the  iterative  inversion  without  assuming  the  off  diagonal  matrix 
elements  are  equal,  (B.17)  (See  section  B  for  more  detail),  can  be  written  as: 

e 

0.5  >  MAX£,e'  tTjnii’m'*/l^emii'm  *  "  0.002*2  ^’(  1/m) 

m=l 

■i 

=0.004/  dx/x  (7.2) 


=  0.004  ) 


0.5  >  0.02 
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where  we  are  using  the  worst  case  conditions  ofA»«.'=  i  and  m=0.  Notice  that 
the  diagonal  dominance  criterion  is  more  than  met.  Ttc  inclination  of  the 
orbit  an  angle  z  from  a  polar  orbit  will  not  change  this  result  since  we  did 
not  assume  any  particular  w{  e)  for  the  infinite  trajectory  length  case  to 
obtain  the  relative  size  of  the  step  function  in  cji. 

Off  diagonal  T-matrix  terms  in  m  will  be  produced  by  higher  order 
spherical  harmonic  terms  in  the  gravity  potential  through  their  effect  on  the 
trajectory  estimate  r  (t).  (In  our  angular  integration  model,  this  will 
produce  <f>  dependence  in  w(TJ  through  rQ(t)*s  influence  on  n(^.)  Terms  in  the 
gravity  field  with  the  index  ±jn'  will  cause  coupling  terms  in  T  between  m  and 
m-an'.  However,  all  gravity  potential  spherical  harmonic  terms  with  !m'!>0  are 
less  than  3x10  “®  of  the  1/r  term  (Gray,  1972)  so  we  expect  the  relative  size 
of  terms  off  diagonal  in  m  due  to  this  effect  to  be  on  the  order  of  3x10  or 
less.  (B.18),  as  mentioned  a  bo  ve^  states  that  the  iterative  inversion  will 
converge  if  the  relative  size  of  the  off  diagonal  elements  is  less  than  0.001. 
Thus,  this  effect  should  not  prevent  the  full  matrix  from  converging 
i terati vely. 

Another  effect  which  must  be  considered  is  roll  and  yaw  of  the 

satellites  with  frequencies  equal  to  nv/(2Tr)  (n=  1,2,3, . )  producing  T 

terms  between  m  and  m+m'.  Because  the  orbit  is  polar,  however,  this  would  not 
produce  terms  coupling  different  m's  but  terms  coupling  different  t's.  Also, 
since  this  effect  would  produce  range  rate  data  which  would  alias  the  effect 
of  gravity  terms,  GSFC  and  APL  have  taken  great  pains  to  ensure  that  this 
effect  will  be  measured  and  taken  out  of  the  range  rate  data  before  further 
analysis. 

Lunar  and  Solar  effects  are  also  too  small  to  effect  the 
convergence  of  the  iterative  inversion  of  the  T-matrix  and  they  will  usually 
be  removed  from  the  range  rate  residual  by  the  least  squares  fit  which 
determines  the  approximate  trajectory. 


8.  CONCLUSIONS 

The  approach  discussed  in  this  paper  of  examining  the  symmetries  of  an 
operator  representation  of  the  observable  proved  very  fruitful  in  simplifying 
the  solution  of  the  matrix  equation  for  determining  the  gravity  potential  from 
low-low  satellite  range  rate  data.  This  still  leaves  the  even  larger  problem 
of  determining  the  matrix  elements  of  the  equation.  This  problem  may  also  be 
simplified  by  utilizing  the  T-matrix  in  the  circular  orbit  infinite  trajectory 
approximation,  which  can  be  very  quickly  calculated  (See  Appendix  A.),  as  a 
zeroth  order  approximation  for  the  actual  matrix.  Further  study  along  these 
lines  will  probably  prove  very  promising.  In  addition,  the  circular  orbit 
infinite  trajectory  approximation  should  prove  very  useful  in  studying  the 
effects  of  various  system  errors  and  interactions  on  the  errors  in  the 
coefficients  of  the  spherical  harmonic  expansion  of  the  gravity  potential. 

This  approach  may  also  nrove  very  fruitful  in  simplifying  other  problems  in 
satellite  geodecy  and  orbit  determination. 
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APPENDIX  A.  CALCULATING  THE  T-MATRIX  IN  THE  CIRCULAR  ORBIT  APPROXIMATION 


The  effect  of  P  on  Yj^^can  be  calculated  by  using  the  following 
formulas  (Messiah,  1966, V.l,  Appendix  B): 

a/ae  =  {l/2)(e'^^L^  -ei^L_)  (A.l) 

where: 

L  =e  +  +  icot(0)  fL  )  (A. 2) 

^  36  3<|) 

and  where: 

L  Y  =  (£(£+1)  -  m(mil))*  Y  (A. 3) 

±  2.  m  2»m±l 

If  we  now  turn  the  full  formalism  of  quantum  mechanics,  we  can  obtain 
a  representation  of  the  T-operator.  In  bra-ket  notation  (Messiah, 1966) ,  one 
can  write  the  matrix  elements  of  T  as: 

=  i  E  <2mlP^w(8)Pl2'm*>  (A. 4) 


where  the  superscript  ^  means  the  hermetlan  adjoint  of  the  operator  and  where 
the  w(e)  is  put  between  P^  and  P  to  indicate  that  the  P  operators  only  operate 
on  the  Yj^^  functions.  Thus  the  T  operator  is: 

T=  iz  (e®®\^--l)w(0)(e-^''^®-l)  (A. 5) 


where  the  hermetian  adjoint  of  3/ae  is  written  as  06\6  to  emphasize  that  it 
operates  in  the  backward  direction. 

If  we  expand  the  exponentials,  collect  terms  of  the  same  order  in 
p,  and  perform  the  sum  over  the  two  directions  of  the  orbit,  we  obtain: 


T  = 

The  lowest  order  term 


«  2n-i 
E  E 
n=l  j=l 

in  p,  Tx 


w(e)(4^)j 
(2n-j)!  j! 
is  given  as: 


=  (  96\6)P^<{e)(  3/30). 


(A. 6) 


(A. 7) 


Setting  w'(r)*l,  we  can  write  w(0)  as: 
w(  0)  =  Nrev/(  iTsi  n0). 


(A. 8) 


From  a  recursion  relation  derived  from  the  fact  that  Y^^^^is  prop0n;ional  to 
e~^^sin9  (Mertzbacher,  1961),  one  can  show  that  the  "operator"  l/sin0 


operating  on  Y  yields: 


where  F 


n 

m 


Y  /sine  »  (22+1)^  i"’  *^2 j 

j=o  t  -n^ — 

^2j+2 

the  number  of  permutations  of  m  object  out  of  n,  is: 

fJJ  =  n(n-l)(n-2) . (n-m+1)  (A. 10) 


■li  Y  , 


(A. 9) 


and  where: 


4  -  2j^  '  1  L  Im+ll  (A. 11) 

Using  (A.l),  (A. 3),  (A. 6),  (A. 9),  and  the  orthonormality  condition 
(6.6),  a  closed  formula  can  be  written  for  the  T-matrix  elements  as  a  sum  of 
coefficients  which  are  calculable  functions  of  a.m.t'.m',  and  p  without  having 
to  perform  any  numerical  integrations.  This  could  be  set  up  in  recursive  form 
on  a  computer  and  accomplished  without  using  a  great  deal  of  computer  time. 
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APPENDIX  8.  SOLVING  THE  MATRIX  EQUATIONS 


B.l  DIRECT  METHODS 


Before  proceeding  it  is  useful  set  up  some  simplifying  termi nolojy.  The 
matrix  equation  to  be  solved  is  given  by: 


TX  =  B 


(B.l) 


where  T  and  B  are  known  and  X  is  to  be  solved  for.  T  is  a  square  matrix  of 
size  NxN,  B  is  an  NxM  matrix  (a  column  vector  when  M=l),  and  the  solution  X  is 
an  NxM  matrix.  In  the  case  of  finding  the  solution  of  the  inversion  of  T,  B  is 
the  NxN  identity  matri)^!,  and  X=T"'.  When  performing  certain  operations  on 
the  equation  in  finding  the  solution,  .t  is  useful  to  introduce  the  concept  of 
the  augmented  matrix,  T^.  T^  is  an  Nx(N+M)  matrix  formed  by  augmenting  the  N 
columns  of  T  with  the  M  columns  of  B;  that  is: 

T^  =  T  (+)  B  (B.2) 


where  (+)  indicates  the  dimensional  addition  of  the  columns  of  T  and  B 
together. 

The  various  direct  numerical  methods  used  to  solve  matrix  equations 
and  invert  matrices  fall  into  3  basic  catagories.  Triangular  Gaussian 
elimination  methods  (Dool  ittle,Crout,  etc.).  Full  Gaussian  elimination 
(Gauss-Jordan),  and  sequential  transformation  methods  (LU  decomposti on, 
partitioned  matrices.  Householder)  (Riess,  1981;  Gastinel,  1970;  Detoma  and 
Wardrip,  1982;  Ralston  and  Rabinowitz,  1978).  The  sequential  transformation 
methods  are  usually  used  only  for  special  cases  to  take  advantage  of  special 
symmetries  such  as  large  blocks  of  zero  elements  (LU  and  partitioning)  or  when 
only  a  small  amount  of  memory  is  available  (Householder).  For  the  general  case 
on  large  computers,  usually  one  of  the  first  two  catagories  is  used. 


Both  of  the  first  two  catagories  rely  on  Gaussian  elimination  (Riess, 
1981;  Gastinel,  1970).  This  consists  of  turning  the  off  diagonal  elements  of 
pivotal  column  (full),  1,  or  a  section  of  the  pivotal  column  (triangular)  of  T 
into  zeros  by  operating  on  T^with: 

=  T^^"^  -  (T?"^  /T^"^  )T^5"^  (8.3) 

U  ij  Ij 


for  all  j  not  equal  ^  in  the  Jordan  method  and  for  c  11  j  greater  than  ^  in  the 
triangular  method.  This  is  performed  N-1  times  witl,  1  =  1  to  N-1  to  produce 
I(+)X  in  the  Jordan  method  and  T"(+)3'  in  the  triangular  method  where  T"  is  an 
upoer  triangular  matrix.  In  the  Jordan  method  the  row  j=l  of  T^is  also 
operated  on  by  dividing  the  1th  row  by  T,j.(The  complication  of  pivoting  to 
minimize  the  error  will  not  be  discussed  here  (Riess,  1981;  Gastinel,  1970).) 
In  the  triangular  method,  X  is  found  by  back-solving  the  modified  equation  as 
follows:  Since  T"  is  upper  triangular,  the  lowest  row  of  the  equation  is  of 
the  '^orm  T"  X  so  i  t  can  be  solved  for  X^ .  by  dividing  Bfj,-  by  Tjuj^  . 

The  value  of  X  •^can  now  be  used  in  the  next  lowes"^  row  of  equation  to  find 


X|\i_i  H  ,  and  this  in  turn  can  be  used  to  find  X.,  ,  ,,  etc.,  until  t^e  whole 
e'quit'ion  is  solved.  ’J 


For  M=l,  and  T  symmetric,  the  number  of  operations  to  perform 
triangular  Gaussian  elimination  (The  triangular  method  takes  fewer  or  as  many 
operations  as  the  Jordan  method  for  all  cases  (Riess,  1981).)  is  approximately 
{1/3)N^  real  floating  point  operations  (flop)  (Riess,  1981;  Gastinel,  1970). 
When  M=N  and  T  is  symmetric  (inversion)  it  takes  approximately  (4/3)N^flop. 

Using  this  method  to  invert  T,  for  T  complex, it  takes  4  times  as  many 
real  operations  as  for  T  real.  Thus  for  complex  and  syitmetric, 

(4/3)N^  (real)  flop  are  required  to  find  A  and  (16/3)N^  (real)  flop  are 
required  to  invert  T. 

It  is  not  enough  just  to  show  that  large  numbers  of  the  matrix 
elements  are  zero  or  near  zero;  when  a  matrix  is  inverted  using  one  of  the 
various  Guassian  elimination  techniques  zeros  and  near  zeros  can  be 
back-filled  with  non-small  values  very  quickly  during  the  factoring  stage 
unless  there  is  some  symmetry  to  prevent  this  from  happening  (Riess,  1981; 
Ouff,  1981).  To  see  how  fast  the  non-zeros  turn  into  zeros,  consider  the 
following.  Suppose  we  have  a  large  WxN  matrix  with  a  relatively  small  number, 
k,  of  non-zero  off  diagonal  elements  in  each  row^with  the  non-zero  elements 
occuring  in  a  random  pattern  from  row  to  row.  In  the  factoring  process  of 
Guassian  elimination  each  row  of  the  matrix  T  other  than  a  pivotal  row,  1,  is 
transformed  by  the  relation  (8.3).  If  T.(^)  were  zero  or  near  zero  and 
and  T|j(n)  were  non-small,  one  can  see  that  T  •  ('^'*'^)would  become  a 
non-small  number.  This  is  what  is  called  back-filling.  After  operating  on  the 
matrix  elements  with  1  pivotal  rows  and  columns,  because  the  probability  of 
T.(^)  being  zero  is  virtually  1  when  there  are  only  a  small  number  of 
n(3rt-zero  off  diagonal  elements,  one  can  show  that  the  average  number  of 
non-zero  terms  in  the  remaining  columns  grows  approximately  as  k2l  until  the 
number  of  non-zero  elements  per  row  in  the  remaining  columns  becomes 
comparable  to  N.  To  keep  this  from  happening,  the  zeros  in  one  row  must  be 
correlated  to  the  zeros  in  another  row  to  a  very  high  order  and  must  remain 
correlated  during  the  elimination  operatati ons. 

3.2  ITERATIVE  METHODS 

The  two  principal  iterative  methods  are  the  Jacobi  method  and  the 
Gauss-Seidel  method  (Riess,  1981;  Gastinel,  1970).  The  Jacobi  method  finds  the 
(n  +  l)th  estimate  of  X  from  the  (r^h  estimate  by; 

=  (8  -  z  T.  ,x('^)  )/T  ,  ,  (B.4) 

1m  im  ij  j  n 

The  Gauss-Seidel  method  uses  the  most  current  estimate  for  which  for  k<i 

is  Each  iteration  takes  2 N"*"  flop  for  M=1  and  for  inverting  a 

matrix  when  T  is  symetric.  Thus,  there  is  no  gain  in  using  tnis  technique  for 
inverting  matrices,  but  there  is  a  very  substantial  gain  for  M=1  as  long  as 
the  solution  converges  in  a  relatively  small  number  of  steps. 

3.3  THE  CO'JVEDGENCE  OF  ITERATIVE  SOLUTIONS 

If  T  is  diagonally  dominant  and  non-singular,  the  Jacobi, 
uau s s - Sei de 1 ,  and  other  similar  iterative  methods  for  finding  the  inverse  of  a 
matrix  or  solving  matrix  equations  will  converge  regardless  of  the  first  guess 


(B.5) 


I  7 


(Riess,  1981;  Gastlnel,  1970).  Diagonal  dominance  is  defined  by: 

I  >Max, 

j 

where  V  means  the  sum  over  j  with  j  not  equal  to  i  and  Max ^  means  the  largest 
value  ^^or  i  =  1  to  N.  In  fact  if  we  define  q  as: 

q  =  Max^  E'IT (B.6) 

j 

one  can  show  that  after  n  iterations: 

llAX^"^  n/l|AX^°^l  <  q"  (B.7) 


where  ||AX^'^^)|  is  the  largest  error  in  the  nth  estimate  of  the  solution,  X. 
This  means  that,  starting  with  as  a  first  guess  for  X,  and  given  q  =1/2, 

the  matrix  elements  would  converge  to  2"'^  times  the  largest  element  in  n 
iterations  or  approximately  to  a  precision  of  2 in  n  iterations.  If  the 
matrices  are  complex,  as  stated  above  each  iteration  in  the  Jacobi  or 
Gauss-Seidel  methods  requires  8N^{real)  flop.  Thus,  obtaining  a  result  with 
2"'^preci sion  requires  only  BN^xn  flop.  For  64  bit  precision,  this  becomes 
512N2  (real)  flop. 

Diagonal  dominance  is  a  stringent  condition  for  large  matrices.  If  all 
the  of  diagonal  elements  were  the  same,  (B.5)  would  become: 


IT^j  I/IT^^  1  =  q/(N-l) 


For  N=41,253  and  q  =1/2,  this  condition  is 


IT.  ./T. .  I 
ij  n 


1.2x10-5. 


(B.3) 


B,4  ITERATIVE  METHOD  FOR  DETERMINING  THE  EXACT  COEFFICIENTS. 

In  this  section,  we  derive  an  iterative  method  which  solves  for  A 
using  X',  the  inverse  of  T' ,  where  T'  is  the  T-matrix  with  all  elements  off 
diagonal  in  m  set  to  zero.  The  method  produces  A  to  64  bit  precision  in  less 
than  1x10^^  (real)  flop.  (See  Sections  B.5  and  7  for  proof  of  convergence.) 
It  can  produce  A  to  any  desired  precision  with  corespondi ngly  fewer  or  more 
flop. 


If  we  observe  that: 

T  =  T'  +  AT, 


(B.9) 


we  can  write: 


TA  =  (T'+  AT)A  =  L. 

Multiplying  both  sides  of  (B.  9)  by  X’,  we  obtain: 


(B.IO) 


(I  +X'  aT)A  -  X’L. 


(B.ll) 


This  equation  can  now  be  solved  by  the  Gauss-Seidel  or  Jacobi  methods  if  the 
matrix  elements  of  T  off  diagonal  in  m  are  small  since  I+X'aT  has  no  matrix 
elements  diagonal  in  m  which  are  off  diagonal  in  i.  (See  Sections  B.5  and  7 
for  proof  of  convergence. ) 

Calculating  A  from  (B.ll)  is  dominated  by  the  number  of  operations  in 
the  multiplication  of  X'  by  aT  and  in  the  iterative  proceedure.  Since  X'  is 
diagonal  in  m,  and  both  X'and  aT  are  symmetric,  it  takes  8N  /3=  5.1x10^^ 
(real)  flop  to  perform  the  multiplication.  Each  iteration  of  the  solution 
takes  8N^  (real)  flop.  Assuming  the  solution  will  converge  64  iterations  or 
less  (See  Section  B.3.),  the  iterative  solution  will  take  SIZN^  =5.5x10^^ 
(real)  flop.  This  means  that  the  total  number  of  operations  to  find  A  is 
1.06x10^^  (real)  flop.  This  is  a  factor  of  88  fewer  flop  than  our  original 
estimate  based  on  (4/3)41,253  ^=9.4x10^^  flop. 

8.5  CONVERGENCE  CRITERION  FOR  MATRIX  ELEMENTS  OFF  DIAGONAL  IN  m 

Based  on  the  iterative  method  of  the  previous  section,  iterative 
convergence  criteria  for  the  T-matrix  elements  off  diagonal  in  m  will  be 
deri ved. 


Since  T  has  no  diagonal  terms  in  m  and  X'  has  only  diagonal  terms  in 
m,  we  note  that  I  +X'  aT  is  equal  to  I  for  diagonal  terms  in  m,and  equal 
to  X*  AT  for  off  diagonal  terms  in  m.  The  criterion  for  convergence  of  the 

iterative  solution  of  (B.ll)  is  that  I+X' aT  be  diagonally  dominant  as  given  by 
(B.5)  and  (B.6),  which  in  this  case,  with  the  noted  properties  of  I+X' 4T, 
become; 


Noting  that: 


1  >  q  =  Max  Z’E'I  £  X.  -  AT,,  ,  , 
tm  2,  m  ^  f mim  2’'m2  m 


(B.12) 


I  |T|  I  =  Max.  z'lT  . 

’  j  iJ 


(B.13) 


is  a  matrix  norm  of  T  which  satisfies  the  Cauchy-Schwarz  inequality  (Riess, 
1982)  ^  w  e  *1  w/  r  I  C  i 

MX'  aTH  <  Mx'M  MaTM  .  (B.14) 

T  li  L/i  ,  CO>i^)exV>o( 

q  <  MX'M  MaTM.  (B.15) 

But  IIX'M  is  on  the  order  of  1/MT'M,  so  we  can  write  the  approximate 
relati  on: 


q  <  MAX 


n,iU 


-1/ 


m 


(B.16) 


Assuming  (T 
we  expect  ti 

can  "divide  out"  the  sum  over  in  the  numerator.  This  yields  the  following 
aporoximate  relation: 


m. 


is  on  the  order  or  bigger  than 
off  diagonal  elements  to  get  smaller  as  m 


m  ' 
and 


(For  this  matrix, 
m'  di ffer. ) .  we 


r 


(B.17) 
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q  <  MAX  .  j;.|T  ,  ,|/!T  , 


where  MAXj^  j^ijpeans  the  maximum  value  ranging  over  allV  and  «•'.  (B.17), 
though  of  ioifiewhat  dubious  derivation,  is  probably  good  for  order  of  magnitude 
calculati ons. 


A  crude  estimate  for  diagonal  dominance  can  be  obtained  from  (B.17)  by 
setting  all  the  off  diagonal  elements  equal.  Setting  q  =0.5,  for  the  worst 
case  of  i  and  equal  to  180  and  m  equal  to  0,  we  obtain: 


•  T IT 1  <  q/(2 

{.mu  m  «,m£  m 


z  ^-l)  =  0.001 
m 


where  the  value  given  is  for  q=  0.5. 


(B.18) 


B.6  PERTURBATION  THEORY  METHOD  FOR  DETERMINING  THE  COVARIANCE  MATRIX 


From  (4.16),  the  covariance  matrix  can  be  determined  if  we  determine 
T"^  and  To  determine  T“^  for  this  purpose,  we  don't  need  an  exact 
solution,  but  only  an  approximate  one  good  enough  for  error  estimates.  The 
following  method  will  obtain  an  approximate  solution  for  T"^  in  only  5x10^^ 
(real )  flop. 

We  want  the  solution  to: 

TX=I  (B.19) 

where  I  is  the  identity  matrix. 

Suppose  we  have  a  solution: 

T'X'=I  (B.20) 

where  T  is  close  to  T' .  Then  by  manipulating  (B.19)  and  (B.20),  we  can 
generate: 


X-X'  =  (I  -  X'T)T~'  (8.21) 

If  we  use  Xj  the  separation  T  =  T'  +  aT  from  Section  B.4^  and  make  the 
approximate  substitution  T  "^  =  X',  we  obtain  the  following  approximate  formula 
(1st  order  perturbation  solution): 

T’'  =  X'  -  X'aTX'  (B.21) 

Given  both  X'aT  and  X'  (from  the  solution  to  TA=L  in  Section  B.4),  it 
takes  8N  2,^  /3  (real)  flop  to  produce  X'TX'  and  4  /3  (real)flop  to 

perform  the  remaining  subtraction.  The  actual  x^^an  oe  calculated  from  (4.17) 
using  y*  defined  by  (4.18)  with  the  integral  in  (4.18)  changed  back  to  the 
integral  over  the  trajectory.  For  3.9x10  ® data  points,  calculating  y^ 
takes  8x10'^  (real)  flop  and  calculating  A^  takes  only  3N=2. 6x10  ^  (real )  flop. 
Thus,  it  takes  approximately  8N  /3  (real)  =  5.1x10  “  (real)  flop  to 
produce  this  estimate  of  the  covariance. 
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REPLY  TO 

ATTN  OF  LWG/Lt  Warner 

SUBJECT  Minutes,  Gravity  Gradiometer  Conference,  12-13  February  1985 


Conference  Attendees 

1.  Please  find  attached  a  copy  of  the  meeting  minutes,  presentations 

(if  individually  requested),  and  list  of  conference  attendees.  The  minutes 
were  compiled  by  Dr.  Chris  Jekeli  and  me  from  abstracts  submitted  by 
various  presenters  and  notes  taken  at  the  conference.  In  some  cases  our 
notes  may  have  been  sketchy  resulting  in  our  interpretation  of  questions 
and  responses  after  each  presentation.  In  all  cases  we  attempted  to 
report  each  presentation  correctly  and  adequately. 

2.  I  have  reserved  the  same  facilities  at  the  USAF  Academy  for  next 
year's  conference  to  be  held  11-12  February  1986.  We'll  start  requesting 
papers  late  this  summer.  Registration  will  proceed  soon  after. 

3.  Thank  you  for  your  interest  in  gravity  gradiometry.  We  look  forward 
to  your  continuing  support. 

DONNA  L.  WARNER  1  Atch:  Minutes 

Conference  Coordinator 


